C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C INCLUDE FILE: mgdraw.inc C C this is a custom include file containing the common blocks shared by the other user routines. C It is not a standard FLUKA routine but has been introduced for programming convenience. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c idcurr = position of the current particle in the particle common block c if idcurr<=0 -> new particle to be followed c integer idbflg, outunit integer tarfrag, fragtrig, stcfrag, bmnfrag, tgtfrag, airfrag integer elsfrag logical trig_stc, trig_bmn, trig_tgt, trig_air, trig_els integer numev integer idcurr integer JSTART, JRUN, JEND parameter (JSTART=1, JRUN=0, JEND=-1) double precision Ethrdep common /comunico/Ethrdep, idbflg, numev, & idcurr, outunit, stcfrag, bmnfrag, tgtfrag, airfrag, & tarfrag, elsfrag, fragtrig, trig_stc, trig_bmn, trig_tgt, & trig_air, trig_els c c particle common block c c nump = number of tracks c c for each of the produced track we register: c c idpa = index in the part common of the track parent c icha = charge c iba = barionic number c idead = with which interaction the particle dies c jpa = kind of particle c igen = kind of interaction c vxi, vyi, vzi = production position of the track c vxf, vyf, vzf = final position of the track c px,py,pz, = production momentum of the track c amass = track mass c tempo = production time of the track c tof = tof of the particle c trlen = track lenght of the particle c integer nump, maxnump parameter(maxnump=8000) integer idpa(maxnump), igen(maxnump) integer icha(maxnump), numreg(maxnump), iba(maxnump) integer idead(maxnump), jpa(maxnump) real vxi(maxnump),vyi(maxnump),vzi(maxnump) real vxf(maxnump),vyf(maxnump),vzf(maxnump) real px(maxnump),py(maxnump),pz(maxnump) real pxf(maxnump),pyf(maxnump),pzf(maxnump) real amass(maxnump), tempo(maxnump), tof(maxnump) real trlen(maxnump) common /particle_common/ vxi, vyi, vzi, & vxf, vyf, vzf,px, py, pz, pxf, pyf, pzf, & amass, tempo, tof, trlen, nump, idpa, igen, & icha, numreg, iba, idead, jpa c ************************************************************* c common needed to score the hit in the start counter c integer nSTC, maxSTC parameter (maxSTC=200) integer idSTC (maxSTC) real xinSTC(maxSTC), xoutSTC(maxSTC), & yinSTC(maxSTC), youtSTC(maxSTC), & zinSTC(maxSTC), zoutSTC(maxSTC) real pxinSTC(maxSTC), pxoutSTC(maxSTC), & pyinSTC(maxSTC), pyoutSTC(maxSTC), & pzinSTC(maxSTC), pzoutSTC(maxSTC) real deSTC(maxSTC), timSTC(maxSTC), & alSTC(maxSTC) common /STC_common/ xinSTC, xoutSTC, yinSTC, & youtSTC, zinSTC, zoutSTC, pxinSTC, pxoutSTC, & pyinSTC, pyoutSTC, pzinSTC, pzoutSTC, & deSTC, alSTC, timSTC, nSTC, idSTC c ************************************************************* c common needed to score the hit in the beam monitor c integer nBMN, maxBMN parameter (maxBMN=1000) integer idBMN (maxBMN), iviewBMN(maxBMN),ilayBMN(maxBMN), & icellBMN(maxBMN) real xinBMN(maxBMN), xoutBMN(maxBMN), yinBMN(maxBMN), & youtBMN(maxBMN), zinBMN(maxBMN), zoutBMN(maxBMN) real pxinBMN(maxBMN), pxoutBMN(maxBMN), pyinBMN(maxBMN), & pyoutBMN(maxBMN), pzinBMN(maxBMN), pzoutBMN(maxBMN) real deBMN(maxBMN), timBMN(maxBMN), alBMN(maxBMN) common /BMN_common/ xinBMN, xoutBMN, yinBMN, youtBMN, zinBMN, & zoutBMN, pxinBMN, pxoutBMN, pyinBMN, pyoutBMN, pzinBMN, & pzoutBMN,deBMN, alBMN, timBMN, nBMN, idBMN, iviewBMN, & ilayBMN, icellBMN c c common for the indexing of the beam monitor c integer MAXBMNREG parameter (MAXBMNREG=200) integer ireg2layBMN(MAXBMNREG), ireg2viewBMN(MAXBMNREG), & ireg2cellBMN(MAXBMNREG) common /BMN_index_common/ireg2layBMN, ireg2viewBMN, ireg2cellBMN c ************************************************************* c common needed to score the hit in the vertex c integer nVTX, maxVTX parameter (maxVTX=300) integer idVTX (maxVTX), ilayVTX(maxVTX) real xinVTX(maxVTX), xoutVTX(maxVTX), yinVTX(maxVTX), & youtVTX(maxVTX), zinVTX(maxVTX), zoutVTX(maxVTX) real pxinVTX(maxVTX), pxoutVTX(maxVTX), pyinVTX(maxVTX), & pyoutVTX(maxVTX), pzinVTX(maxVTX), pzoutVTX(maxVTX) real deVTX(maxVTX), timVTX(maxVTX), alVTX(maxVTX) common /VTX_common/ xinVTX, xoutVTX, yinVTX, youtVTX, & zinVTX, zoutVTX, pxinVTX, pxoutVTX, pyinVTX, pyoutVTX, & pzinVTX, pzoutVTX, deVTX, alVTX, timVTX, nVTX, idVTX, & ilayVTX c ************************************************************* c common needed to score the hit in the inner tracker c integer nITR, maxITR parameter (maxITR=300) integer idITR (maxITR), isensITR(maxITR) c integer idITR (maxITR), iplumeITR(maxITR), ilayITR(maxITR), c & imimoITR(maxITR) real xinITR(maxITR), xoutITR(maxITR), yinITR(maxITR), & youtITR(maxITR), zinITR(maxITR), zoutITR(maxITR) real pxinITR(maxITR), pxoutITR(maxITR), pyinITR(maxITR), & pyoutITR(maxITR), pzinITR(maxITR), pzoutITR(maxITR) real deITR(maxITR), timITR(maxITR), alITR(maxITR) common /ITR_common/ xinITR, xoutITR, yinITR, youtITR, & zinITR, zoutITR, pxinITR, pxoutITR, pyinITR, pyoutITR, & pzinITR, pzoutITR, deITR, alITR, timITR, nITR, idITR, & isensITR c & iplumeITR, ilayITR, imimoITR c c common for the indexing of the ITR c integer MAXITRREG parameter (MAXITRREG=300) integer ireg2sensITR(MAXITRREG) common /ITR_index_common/ireg2sensITR c integer ireg2plumeITR(MAXITRREG), ireg2layITR(MAXITRREG), c & ireg2mimoITR(MAXITRREG) c common /ITR_index_common/ireg2plumeITR, ireg2layITR, c & ireg2mimoITR c ************************************************************* c common needed to score the hit in the microstrip detector c integer nMSD, maxMSD parameter (maxMSD=1000) integer idMSD (maxMSD), ilayMSD(maxMSD) real xinMSD(maxMSD), xoutMSD(maxMSD), yinMSD(maxMSD), & youtMSD(maxMSD), zinMSD(maxMSD), zoutMSD(maxMSD) real pxinMSD(maxMSD), pxoutMSD(maxMSD),pyinMSD(maxMSD), & pyoutMSD(maxMSD), pzinMSD(maxMSD), pzoutMSD(maxMSD) real deMSD(maxMSD), timMSD(maxMSD), alMSD(maxMSD) common /MSD_common/ xinMSD, xoutMSD, yinMSD, youtMSD, zinMSD, & zoutMSD, pxinMSD, pxoutMSD, pyinMSD, pyoutMSD, pzinMSD, & pzoutMSD,deMSD, alMSD, timMSD, nMSD, idMSD, & ilayMSD c c common for the indexing of the microstrip detector c integer MAXMSDREG parameter (MAXMSDREG=300) integer ireg2layMSD(MAXMSDREG) common /MSD_index_common/ireg2layMSD c c ************************************************************* c common needed to score the hit in the scintillator c integer nSCN, maxSCN parameter (maxSCN=1000) integer idSCN(maxSCN) integer istripSCN(maxSCN), iviewSCN(maxSCN) real xinSCN(maxSCN), xoutSCN(maxSCN), yinSCN(maxSCN), & youtSCN(maxSCN), zinSCN(maxSCN), zoutSCN(maxSCN) real pxinSCN(maxSCN), pxoutSCN(maxSCN), & pyinSCN(maxSCN), pyoutSCN(maxSCN), & pzinSCN(maxSCN), pzoutSCN(maxSCN) real deSCN(maxSCN),alSCN(maxSCN), timSCN(maxSCN) common /SCN_common/ xinSCN, & xoutSCN, yinSCN, youtSCN, zinSCN, zoutSCN, & pxinSCN, pxoutSCN, pyinSCN, pyoutSCN, pzinSCN, & pzoutSCN, deSCN, alSCN, timSCN, nSCN, idSCN, & istripSCN, iviewSCN c c common for the indexing of the SCN c integer MAXSCNREG parameter (MAXSCNREG=300) integer ireg2stripSCN(MAXSCNREG),ireg2viewSCN(MAXSCNREG) common /SCN_index_common/ireg2stripSCN,ireg2viewSCN c ************************************************************* c common neede to score the hit in the Calo regions c integer nCAL, maxCAL parameter (maxCAL=8000) integer idCAL(maxCAL),icryCAL(maxCAL) real xinCAL(maxCAL), xoutCAL(maxCAL), yinCAL(maxCAL), & youtCAL(maxCAL), zinCAL(maxCAL), zoutCAL(maxCAL) real pxinCAL(maxCAL), pxoutCAL(maxCAL), pyinCAL(maxCAL), & pyoutCAL(maxCAL), pzinCAL(maxCAL), pzoutCAL(maxCAL) real deCAL(maxCAL), timCAL(maxCAL), alCAL(maxCAL) common /CAL_common/ xinCAL, xoutCAL, yinCAL, youtCAL, zinCAL, & zoutCAL, pxinCAL, pxoutCAL, pyinCAL, pyoutCAL, pzinCAL, & pzoutCAL, deCAL, alCAL, timCAL, nCAL, idCAL, icryCAL c c common for the indexing of the SCN c integer MAXCALREG parameter (MAXCALREG=700) integer ireg2cryCAL(MAXCALREG) common /CAL_index_common/ireg2cryCAL c ************************************************************* c common needed to score the boundary crossing info c and put the right pointer to the particle common of the produced electrons c integer nCROSS, maxCROSS parameter (maxCROSS=25000) integer idCROSS(maxCROSS), nregCROSS(maxCROSS) integer nregoldCROSS(maxCROSS) real pxCROSS(maxCROSS), pyCROSS(maxCROSS), pzCROSS(maxCROSS) real xCROSS(maxCROSS), yCROSS(maxCROSS), zCROSS(maxCROSS) real tCROSS(maxCROSS), amaCROSS(maxCROSS), chCROSS(maxCROSS) common /crossing_common/ & pxCROSS, pyCROSS, pzCROSS, & xCROSS, yCROSS, zCROSS, tCROSS, chCROSS, amaCROSS, nCROSS, & idCROSS, nregCROSS, nregoldCROSS c ************************************************************ c Numbers of fluka region to be scored c integer nregaria, nregMagAir, nregSTC, nregtarg, & nregFirstBMN, nregLastBMN, nregFirstVTX, nregLastVTX, & nregFirstITR, nregLastITR, nregFirstMSD, nregLastMSD, & nregFirstSCN, nregLastSCN, nregFirstCAL, nregLastCAL common/region_common/ & nregaria, nregMagAir, nregSTC, nregtarg, nregFirstBMN, & nregLastBMN, nregFirstVTX, nregLastVTX, nregFirstITR, & nregLastITR, nregFirstMSD, nregLastMSD, nregFirstSCN, & nregLastSCN, nregFirstCAL, nregLastCAL c *********************************************************** c absolute light efficiency of material in quenched formula c double precision abs_STC,abs_BMN,abs_VTX,abs_ITR, & abs_MSD, abs_SCN, abs_CAL parameter (abs_STC=1.0) parameter (abs_BMN=1.0) parameter (abs_VTX=1.0) parameter (abs_ITR=1.0) parameter (abs_MSD=1.0) parameter (abs_SCN=1.0) parameter (abs_CAL=1.0) c *********************************************************** c interaction stack: common block for the parent particle c recording with the info on the interaction point c integer idfluka(maxnump) integer numint, maxint parameter (MAXINT=30000) double precision xint(maxint),yint(maxint),zint(maxint) integer intpa(maxint),intcode(maxint) common /particle_stack/xint,yint,zint,intpa,intcode, & idfluka,numint C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C INCLUDE FILE: parameters .inc C C This is another custom include file containing detectors parameters needed and used by the other user routines. C Since these parameters are strictly related to the geometry, this file is automatically generated when producing C the FLUKA geometry with the MakeGeo macro to reduce the probability of introducing errors C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c BEAM MONITOR PARAMETERS integer ncellBMN parameter (ncellBMN = 3) integer nlayBMN parameter (nlayBMN = 6) c VERTEX PARAMETERS integer nlayVTX parameter (nlayVTX = 4) c INNER TRACKER PARAMETERS integer nsensITR parameter (nsensITR = 32) c MSD PARAMETERS integer nlayMSD parameter (nlayMSD = 6) c MAGNETS PARAMETERS double precision MagCenterX, MagCenterY, MagCenterZ parameter (MagCenterX=0.00000D+00) parameter (MagCenterY=0.00000D+00) parameter (MagCenterZ=19.00000D+00) character*50 mapname parameter (mapname='./data/MagneticMap.txt') c SCINTILLATOR PARAMETERS integer nstripSCN parameter(nstripSCN = 20) c CALORIMETER PARAMETERS integer ncryCAL parameter(ncryCAL = 333) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: usrini.f C C This routine performs the initialization phase. It is called at the beginning of the run by the USRICALL card in the input file. C The meaning of the numerical parameters in this data card is user defined. C Therefore the calling parameters can provide useful paramters or flags to drive the routines actions. C In the FOOT input, the USRICAL card provides two flags: a debug flag producing, if activated, a verbose output, C and an event display flag. Since the customized scoring in the detectors is performed on a region basis, the main task C of the USRINI subroutine is to recognize and store the region names. C In fact, each region of the experimental setup has a double identifier: the first one is its name defined and used by the user, C while the second one is a sequential number internally assigned and used by FLUKA. Hence, an algorithm capable C of linking the region name to the FLUKA internal number has been implemented in this subroutine. C This allow the user to easily recall in the other user routines the regions in which the scoring must be performed. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE USRINI ( WHAT, SDUM ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1991-2005 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * USeR INItialization: this routine is called every time the * * USRICALL card is found in the input stream * * * * * * Created on 01 january 1991 by Alfredo Ferrari & Paola Sala * * Infn - Milan * * * * Last change on 20-mar-05 by Alfredo Ferrari * * * * * *----------------------------------------------------------------------* * DIMENSION WHAT (6) CHARACTER SDUM*8 c INCLUDE '(FLKMAT)' include "../mgdraw.inc" include "../parameters.inc" character*8 REGNAM integer cellBMN, strip, cry, layMSD integer sensITR c integer plmITR, layITR, mimoITR integer ia,ib,ic,id,ie,ig,ih,il,im * * * Don't change the following line: LUSRIN = .TRUE. * *** Write from here on *** * c idbflg = int(what(1)) fragtrig = int(what(2)) Ethrdep = what(3)/1000. c c init of geometry parameters c write(*,*)" " write(*,*)"fragtrig= ",fragtrig," Ethrdep= ",Ethrdep c outunit=60 open(unit=outunit,file='TXT.dat',form='formatted',status='new') c c writing the header c if(idbflg.ne.10) write(outunit,*) int(abs(fragtrig)), Ethrdep c c find the region number of the region of interest c nregtarg = 0 nregaria = 0 nregMagAir = 0 nregSTC = 0 nregFirstVTX = 1000000 nregLastVTX = 0 nregFirstBMN = 1000000 nregLastBMN = 0 nregFirstITR = 1000000 nregLastITR = 0 nregFirstMSD = 1000000 nregLastMSD = 0 nregFirstSCN = 1000000 nregLastSCN = 0 nregFirstCAL = 1000000 nregLastCAL = 0 c ia = 1 ib = 1 ic = 1 id = 1 ie = 1 ig = 1 ih = 1 ik = 1 il = 1 im = 1 c do nn = 1,MAXBMNREG ireg2viewBMN(nn) = -10 ireg2cellBMN(nn) = -10 ireg2layBMN(nn) = -10 end do do nn = 1,MAXITRREG ireg2sensITR(nn) = -10 end do do nn = 1,MAXMSDREG ireg2layMSD(nn) = -10 end do do nn = 1,MAXSCNREG ireg2stripSCN(nn) = -10 ireg2viewSCN(nn) = -10 end do do nn = 1,MAXCALREG ireg2cryCAL(nn) = -10 end do c c do ii = 1,NREGS call GEOR2N ( ii, REGNAM, IERR ) if(ierr.eq.0) then if(REGNAM.eq.'AIR1') then nregaria=ii else if(REGNAM.eq.'MAG_AIR') then nregMagAir=ii elseif(REGNAM.eq.'STC')then nregSTC=ii elseif(REGNAM(1:5).eq.'BMN_C') then read(REGNAM(7:8),*) cellBMN if (REGNAM(6:6).eq.'0')then ireg2viewBMN(ii) = 1 ireg2layBMN(ii) = cellBMN/ncellBMN ireg2cellBMN(ii) = cellBMN - cellBMN/ncellBMN*ncellBMN elseif(REGNAM(6:6).eq.'1') then ireg2viewBMN(ii) = -1 ireg2layBMN(ii) = cellBMN/ncellBMN ireg2cellBMN(ii) = cellBMN - cellBMN/ncellBMN*ncellBMN endif if(ig.eq.1) then nregFirstBMN=ii elseif(ig.eq.(nlayBMN*ncellBMN*2)) then nregLastBMN=ii endif ig = ig + 1 elseif(REGNAM.eq.'TARGET')then nregtarg = ii elseif(REGNAM(1:4).eq.'VTXE') then if(ia.eq.1) then nregFirstVTX=ii elseif(ia.eq.nlayVTX) then nregLastVTX=ii endif ia = ia + 1 elseif(REGNAM(1:4).eq.'ITRE') then read(REGNAM(5:6),*) sensITR ireg2sensITR(ii) = sensITR if(ib.eq.1) then nregFirstITR=ii elseif(ib.eq.nsensITR) then nregLastITR=ii endif ib = ib + 1 elseif(REGNAM(1:4).eq.'MSDS') then read(REGNAM(5:5),*) layMSD ireg2layMSD(ii) = layMSD if(il.eq.1) then nregFirstMSD=ii elseif(il.eq.(nlayMSD)) then nregLastMSD=ii endif il = il + 1 elseif(REGNAM(1:3).eq.'SCN') then if (REGNAM(4:4).eq.'0')then ireg2viewSCN(ii) = 0 elseif(REGNAM(4:4).eq.'1') then ireg2viewSCN(ii) = 1 endif read(REGNAM(5:6),*) strip ireg2stripSCN(ii) = strip if(ic.eq.1) then nregFirstSCN=ii endif nregLastSCN=ii ic = ic + 1 elseif(REGNAM(1:3).eq.'CAL') then read(REGNAM(4:6),*) cry ireg2cryCAL(ii) = cry if(id.eq.1) then nregFirstCAL=ii endif nregLastCAL = ii id = id + 1 endif endif end do c c c write(*,*)'======================================' write(*,*)'USRINI: idbflg = ',idbflg write(*,*)'USRINI: = ',ib, " ", nsensITR c if(((nregtarg*nregLastVTX*nregaria*nregLastSCN*nregSTC*nregLastITR & *nregLastCAL*nregLastBMN & *nregLastMSD).eq.0).or.(nregFirstVTX.eq.1000000).or. & (nregFirstCAL.eq.1000000).or.(nregFirstSCN.eq.1000000).or. & (nregFirstBMN.eq.1000000).or. & (nregFirstMSD.eq.1000000).or. & (nregFirstITR.eq.1000000)) then write(*,*)'I could not find all the regions!!!!' endif write(*,*)'**************** Inizio Geometria *****************' write(*,*)' nregaria = ',nregaria write(*,*)' nregSTC = ',nregSTC write(*,*)' nregFirstBMN = ',nregFirstBMN write(*,*)' nregLastBMN = ',nregLastBMN write(*,*)' nregtarg = ',nregtarg write(*,*)' nregFirstVTX = ',nregFirstVTX write(*,*)' nregLastVTX = ',nregLastVTX write(*,*)' nregFirstITR = ',nregFirstITR write(*,*)' nregLastITR = ',nregLastITR write(*,*)' nregFirstMSD = ',nregFirstMSD write(*,*)' nregLastMSD = ',nregLastMSD write(*,*)' nregFirstSCN = ',nregFirstSCN write(*,*)' nregLastSCN = ',nregLastSCN write(*,*)' nregFirstCAL = ',nregFirstCAL write(*,*)' nregLastCAL = ',nregLastCAL write(*,*)'**************** Fine Geometria *******************' write(*,*)'' write(*,*)' ' write(*,*)'======================================' RETURN *=== End of subroutine Usrini ================================* END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: usrein.f C C It is called at the beginning of each event before the sampled primary is transported. C In FOOT simulations, it allows the user to initialize the event by zeroing the output arrays of the common C blocks declared in the mgdraw.inc file. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE USREIN INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(CASLIM)' * *----------------------------------------------------------------------* * * * Copyright (C) 1991-2005 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * USeR Event INitialization: this routine is called before the * * showering of an event is started, but after the source particles * * of that event have been already loaded on the stack * * * * Created on 01 january 1991 by Alfredo Ferrari & Paola Sala * * Infn - Milan * * * * Last change on 09-apr-99 by Alfredo Ferrari * * * * * *----------------------------------------------------------------------* * c include 'mgdraw.inc' integer ii c if(idbflg .gt. 0)then write(*,*)'***************************************************' write(*,*)' EVENT',ncase write(*,*)'***************************************************' write(*,*)'' endif c numev = ncase tarfrag = 0 c do ii = 1,maxnump if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing nump' idpa(ii) = 0 igen(ii) = 0 icha(ii) = 0 numreg(ii) = 0 iba(ii) = 0 idead(ii) = 0 jpa(ii) = 0 vxi(ii) = 0. vyi(ii) = 0. vzi(ii) = 0. vxf(ii) = 0. vyf(ii) = 0. vzf(ii) = 0. px(ii) = 0. py(ii) = 0. pz(ii) = 0. pxf(ii) = 0. pyf(ii) = 0. pzf(ii) = 0. amass(ii) = 0. tempo(ii) = 0. tof(ii) = 0. trlen(ii) = 0. c idfluka(ii) = 0 ! aux variables for particle latching c end do nump = 0 c c c do ii = 1,maxint if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing numint' xint(ii) = 0 yint(ii) = 0 zint(ii) = 0 intpa(ii) = 0 intcode(ii) = 0 end do numint = 0 c c start counter c nSTC = 0 do ii = 1,maxSTC if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing stc' idSTC(ii) = 0 xinSTC(ii) = 0. xoutSTC(ii) = 0. yinSTC(ii) = 0. youtSTC(ii) = 0. zinSTC(ii) = 0. zoutSTC(ii) = 0. pxinSTC(ii) = 0. pxoutSTC(ii) = 0. pyinSTC(ii) = 0. pyoutSTC(ii) = 0. pzinSTC(ii) = 0. pzoutSTC(ii) = 0. deSTC(ii) = 0. alSTC(ii) = 0. timSTC(ii) = 0. end do c c beam monitor c nBMN = 0 do ii = 1,maxBMN if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing bm' idBMN(ii) = 0 iviewBMN(ii) = 0 ilayBMN(ii) = 0 icellBMN(ii) = 0 xinBMN(ii) = 0. xoutBMN(ii) = 0. yinBMN(ii) = 0. youtBMN(ii) = 0. zinBMN(ii) = 0. zoutBMN(ii) = 0. pxinBMN(ii) = 0. pxoutBMN(ii) = 0. pyinBMN(ii) = 0. pyoutBMN(ii) = 0. pzinBMN(ii) = 0. pzoutBMN(ii) = 0. deBMN(ii) = 0. alBMN(ii) = 0. timBMN(ii) = 0. end do c c vertex c nVTX = 0 do ii = 1,maxVTX if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing vtx' idVTX(ii) = 0 ilayVTX(ii) = 0 xinVTX(ii) = 0. xoutVTX(ii) = 0. yinVTX(ii) = 0. youtVTX(ii) = 0. zinVTX(ii) = 0. zoutVTX(ii) = 0. pxinVTX(ii) = 0. pxoutVTX(ii) = 0. pyinVTX(ii) = 0. pyoutVTX(ii) = 0. pzinVTX(ii) = 0. pzoutVTX(ii) = 0. deVTX(ii) = 0. alVTX(ii) = 0. timVTX(ii) = 0. end do c c inner tracker c nITR = 0 do ii = 1,maxITR if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing itr' idITR(ii) = 0 isensITR(ii)= 0 xinITR(ii) = 0. xoutITR(ii) = 0. yinITR(ii) = 0. youtITR(ii) = 0. zinITR(ii) = 0. zoutITR(ii) = 0. pxinITR(ii) = 0. pxoutITR(ii) = 0. pyinITR(ii) = 0. pyoutITR(ii) = 0. pzinITR(ii) = 0. pzoutITR(ii) = 0. deITR(ii) = 0. alITR(ii) = 0. timITR(ii) = 0. end do c c second drift chamber c nMSD = 0 do ii = 1,maxMSD if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing msd' idMSD(ii) = 0 ilayMSD(ii) = 0 xinMSD(ii) = 0. xoutMSD(ii) = 0. yinMSD(ii) = 0. youtMSD(ii) = 0. zinMSD(ii) = 0. zoutMSD(ii) = 0. pxinMSD(ii) = 0. pxoutMSD(ii) = 0. pyinMSD(ii) = 0. pyoutMSD(ii) = 0. pzinMSD(ii) = 0. pzoutMSD(ii) = 0. deMSD(ii) = 0. alMSD(ii) = 0. timMSD(ii) = 0. end do c c scintillator c nSCN = 0 do ii= 1, maxSCN if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing scn' idSCN(ii) = 0 istripSCN(ii) = 0 iviewSCN(ii) = 0 xinSCN(ii) = 0. yinSCN(ii) = 0. zinSCN(ii) = 0. pxinSCN(ii) = 0. pyinSCN(ii) = 0. pzinSCN(ii) = 0. timSCN(ii) = 0. xoutSCN(ii) = 0. youtSCN(ii) = 0. zoutSCN(ii) = 0. pxoutSCN(ii) = 0. pyoutSCN(ii) = 0. pzoutSCN(ii) = 0. deSCN(ii) = 0. alSCN(ii) = 0. enddo c c calorimeter c nCAL = 0 do ii= 1, maxCAL if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing cal' idCAL(ii) = 0 icryCAL(ii) = 0 xinCAL(ii) = 0. yinCAL(ii) = 0. zinCAL(ii) = 0. pxinCAL(ii) = 0. pyinCAL(ii) = 0. pzinCAL(ii) = 0. timCAL(ii) = 0. xoutCAL(ii) = 0. youtCAL(ii) = 0. zoutCAL(ii) = 0. pxoutCAL(ii) = 0. pyoutCAL(ii) = 0. pzoutCAL(ii) = 0. deCAL(ii) = 0. alCAL(ii) = 0. enddo c c crossings c nCROSS = 0 do ii= 1, nmaxCROSS if (ncase.eq.1 .AND. ii.eq.1) write(*,*) 'zeroing cross' idCROSS(ii) = 0 nregCROSS(ii) = 0 nregoldCROSS(ii) = 0 pxCROSS(ii) = 0. pyCROSS(ii) = 0. pzCROSS(ii) = 0. xCROSS(ii) = 0. yCROSS(ii) = 0. zCROSS(ii) = 0. tCROSS(ii) = 0. chCROSS(ii) = 0. amaCROSS(ii) = 0. end do c c triggers c stcfrag = 0 bmnfrag = 0 tgtfrag = 0 airfrag = 0 elsfrag = 0 trig_tgt = .false. trig_bmn = .false. trig_stc = .false. trig_air = .false. trig_els = .false. c c RETURN *=== End of subroutine Usrein =========================================* END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: mgdraw.f C C This subroutine is the core of the output building, since it handles the energy depositions recording. C It allows to intercept the transport and the interaction processes at every step. This subroutine is activated by C the USRDUMP data card, and it is used to write the output file where all or selected transport events are recorded. C Several entry points can be found in this subroutine: C - MGDRAW, in which the customized scoring in the detector regions is coded and the energy deposition along each step is C calculated; C - SODRAW, which manages the injection of the primary particlein each event; C - BXDRAW, in which the scoring at region crossings is performed; C - ENDRAW, which handles local depositions of particles below threshold; C - USDRAW, where the simulation searches inelastic nuclear interactions in the target or in other regions of interest C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE MGDRAW ( ICODE, MREG ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2006 by Alfredo Ferrari * * All Rights Reserved. * * * * * * MaGnetic field trajectory DRAWing: actually this entry manages * * all trajectory dumping for * * drawing * * * * Created on 01 march 1990 by Alfredo Ferrari * * INFN - Milan * * Last change 05-may-06 by Alfredo Ferrari * * INFN - Milan * * * *----------------------------------------------------------------------* * INCLUDE '(CASLIM)' INCLUDE '(COMPUT)' INCLUDE '(SOURCM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(GENSTK)' INCLUDE '(MGDDCM)' INCLUDE '(PAPROP)' INCLUDE '(QUEMGD)' INCLUDE '(SUMCOU)' INCLUDE '(TRACKR)' INCLUDE '(LTCLCM)' c include 'mgdraw.inc' * DIMENSION DTQUEN ( MXTRCK, MAXQMG ) CHARACTER*20 FILNAM LOGICAL LFCOPE SAVE LFCOPE DATA LFCOPE / .FALSE. / c integer ICPART, IBPART character*8 NEWREGNAM,MREGNAM * CHARACTER*8 MLATNAM logical ldead double precision erawSTC, equenchedSTC double precision erawBMN, equenchedBMN double precision erawMSD, equenchedMSD double precision erawVTX, equenchedVTX double precision erawITR, equenchedITR double precision erawSCN, equenchedSCN double precision erawCAL, equenchedCAL * *----------------------------------------------------------------------* * * * Icode = 1: call from Kaskad * * Icode = 2: call from Emfsco * * Icode = 3: call from Kasneu * * Icode = 4: call from Kashea * * Icode = 5: call from Kasoph * * * *----------------------------------------------------------------------* if(idbflg.lt.0) then IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', & FORM = 'UNFORMATTED' ) END IF WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK), & SNGL (WTRACK) WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)), & SNGL (ZTRACK (I)), I = 0, NTRACK ), & ( SNGL (DTRACK (I)), I = 1, MTRACK ), & SNGL (CTRACK) endif if(idbflg.gt.1) then call GEOR2N ( mreg, MREGNAM, IERR ) write(*,*)' ' write(*,*)'------------- mgdraw: Ev =',ncase,' -------------' write(*,*)'jtrack = ',jtrack,' regione = ',MREGNAM write(*,*)' ' endif c c************* PARTICLES AND HEAVY IONS WITH JTRACK=-2 IF((JTRACK.GE.1).OR.((JTRACK.LE.-2).AND.(JTRACK.GE.-6)))THEN ICPART = ICHRGE(JTRACK) IBPART = IBARCH(JTRACK) AMPART = AM(JTRACK) c*************HEAVY IONS WITH JTRACK.LT.-6 ELSEIF(JTRACK.LT.-6)THEN ICPART = ICHEAV(-JTRACK) IBPART = IBHEAV(-JTRACK) AMPART = AMNHEA(-JTRACK) ELSE ICPART=int(zfrttk) C IBPART=0 C ICPART = ICHRGE(JTRACK) IBPART = IBARCH(JTRACK) AMPART = AM(JTRACK) ENDIF c c************************************************************** CALL GEOR2N(MREG,MREGNAM,IERR) CALL GEON2R(MREGNAM,MREG,IERR) if(idbflg.gt.1) then write(*,*)'ptrack = ',ptrack,' ekin = ',etrack-ampart write(*,*)'zfrttk = ',zfrttk,' icpart = ',icpart,' m =',ampart WRITE (*,*)'x,y,z = ', ( SNGL (XTRACK (I)), SNGL (YTRACK (I)), & SNGL (ZTRACK (I)), I = 0, NTRACK ) WRITE (*,*)'cx,cy,cz = ',SNGL(CXTRCK),SNGL(CYTRCK),SNGL(CZTRCK) WRITE (*,*)'px,py,pz = ',SNGL(ptrack*CXTRCK), & SNGL(ptrack*CYTRCK),SNGL(ptrack*CZTRCK) write(*,*)'time = ',atrack,' step = ',ctrack,' path = ',cmtrck endif call UpdateCurrPart(mreg,icpart,ibpart,ampart,icode, & xtrack(0),ytrack(0),ztrack(0),JRUN) if (idbflg.lt.0) then IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)) & , I = 1, MTRACK ), & JBK = 1, NQEMGD ) ENDIF endif c ***************************************************************** c inside the start counter c if( mreg.eq.nregSTC )then erawSTC = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawSTC = erawSTC + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) do ii = 1,mtrack equenchedSTC = equenchedSTC + dtquen(ii,3) end do equenchedSTC = equenchedSTC*abs_STC endif endif if(erawSTC.gt.0) then call score_STC(mreg,erawSTC,equenchedSTC, & xtrack(0),ytrack(0),ztrack(0),xtrack(ntrack),ytrack(ntrack), & ztrack(ntrack)) endif endif c ***************************************************************** c inside the beam monitor c if( mreg.ge.nregFirstBMN.and.mreg.le.nregLastBMN )then erawBMN = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawBMN = erawBMN + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) do ii = 1,mtrack equenchedBMN = equenchedBMN + dtquen(ii,3) end do equenchedBMN = equenchedBMN*abs_BMN endif endif if(erawBMN.gt.0) then call score_BMN(mreg,erawBMN,equenchedBMN, & xtrack(0),ytrack(0),ztrack(0),xtrack(ntrack),ytrack(ntrack), & ztrack(ntrack)) endif endif c ************************************************************* c inside the vertex c if( mreg.le.nregLastVTX .and. mreg.ge.nregFirstVTX )then erawVTX = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawVTX = erawVTX + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) do ii = 1,mtrack equenchedVTX = equenchedVTX + dtquen(ii,3) end do equenchedVTX = equenchedVTX*abs_VTX endif endif if(erawVTX.gt.0) then call score_VTX(mreg,erawVTX,equenchedVTX, & xtrack(0),ytrack(0),ztrack(0),xtrack(ntrack),ytrack(ntrack), & ztrack(ntrack)) endif endif c ************************************************************* c inside the inner tracker c if( mreg.le.nregLastITR .and. mreg.ge.nregFirstITR )then erawITR = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawITR = erawITR + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) do ii = 1,mtrack equenchedITR = equenchedITR + dtquen(ii,3) end do equenchedITR = equenchedITR*abs_ITR endif endif if(erawITR.gt.0) then call score_ITR(mreg,erawITR,equenchedITR, & xtrack(0),ytrack(0),ztrack(0),xtrack(ntrack),ytrack(ntrack), & ztrack(ntrack)) endif endif c ***************************************************************** c inside the second drift chamber c if( mreg.ge.nregFirstMSD.and.mreg.le.nregLastMSD )then erawMSD = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawMSD = erawMSD + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) c c DTQUEN(MTRACK,1) e' il rilascio di energia quenchato nella camera a drift c do ii = 1,mtrack equenchedMSD = equenchedMSD + dtquen(ii,3) end do equenchedMSD = equenchedMSD*abs_MSD endif endif if(erawMSD.gt.0) then call score_MSD(mreg,erawMSD,equenchedMSD, & xtrack(0),ytrack(0),ztrack(0),xtrack(ntrack),ytrack(ntrack), & ztrack(ntrack)) endif endif c c ********************************************************** c Inside the scintillator c if( mreg.le.nregLastSCN .and. mreg.ge.nregFirstSCN ) then erawSCN = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawSCN = erawSCN + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) do ii = 1,mtrack equenchedSCN = equenchedSCN + dtquen(ii,3) end do equenchedSCN = equenchedSCN*abs_plastic endif END IF if(erawSCN.gt.0) then call score_SCN(mreg,erawSCN,equenchedSCN,xtrack(0), & ytrack(0), & ztrack(0),xtrack(ntrack),ytrack(ntrack),ztrack(ntrack)) endif endif c ************************************************************* c inside the calorimeter c if( mreg.le.nregLastCAL .and. mreg.ge.nregFirstCAL )then erawCAL = 0. IF ( MTRACK .GT. 0 )THEN do ii = 1,MTRACK erawCAL = erawCAL + dtrack(ii) end do IF ( LQEMGD )THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) do ii = 1,mtrack equenchedCAL = equenchedCAL + dtquen(ii,3) end do equenchedCAL = equenchedCAL*abs_CAL endif endif if(erawCAL.gt.0) then call score_CAL(mreg,erawCAL,equenchedCAL,xtrack(0), & ytrack(0), ztrack(0),xtrack(ntrack),ytrack(ntrack), & ztrack(ntrack)) endif endif RETURN * *================================================================* * * * Boundary-(X)crossing DRAWing: * * * * Icode = 1x: call from Kaskad * * 19: boundary crossing * * Icode = 2x: call from Emfsco * * 29: boundary crossing * * Icode = 3x: call from Kasneu * * 39: boundary crossing * * Icode = 4x: call from Kashea * * 49: boundary crossing * * Icode = 5x: call from Kasoph * * 59: boundary crossing * * * *======================================================================* * * ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO ) c************* PARTICLES AND HEAVY IONS WITH JTRACK=-2 IF((JTRACK.GE.1).OR.((JTRACK.LE.-2).AND.(JTRACK.GE.-6)))THEN ICPART = ICHRGE(JTRACK) IBPART = IBARCH(JTRACK) AMPART = AM(JTRACK) c************* HEAVY IONS WITH JTRACK.LT.-6 ELSEIF(JTRACK.LT.-6)THEN ICPART = ICHEAV(-JTRACK) IBPART = IBHEAV(-JTRACK) AMPART = AMNHEA(-JTRACK) ELSE ICPART= ICHRGE(JTRACK) ICPART = ICHRGE(JTRACK) IBPART = IBARCH(JTRACK) AMPART = AM(JTRACK) ENDIF c c debug printing c if(idbflg.gt.0) then call GEOR2N ( newreg, NEWREGNAM, IERR ) call GEOR2N ( mreg, MREGNAM, IERR ) write(*,*)' ' write(*,*)'------------- bxdraw: Ev =',ncase,' -------------' write(*,*)'jtrack = ',jtrack,' icode = ',icode write(*,*)'da -> ',mregnam,'a -> = ',newregnam write(*,*)'zfrttk = ',zfrttk,' icpart = ',icpart, & ' m =',AMPART write(*,*)'ptrack = ',ptrack,' ekin = ',etrack-AMPART write(*,*)'x,y,z = ',XSCO,YSCO,ZSCO WRITE (*,*)'cx,cy,cz = ',SNGL(CXTRCK), SNGL(CYTRCK) & ,SNGL(CZTRCK) write(*,*)'number of track segments: ntrack = ',NTRACK do ii =1,ntrack write(*,*)'track seg # ',ii,' track seg lenght = ' & ,ttrack(ii) end do write(*,*)'number of en deposition: mtrack = ',MTRACK do ii =1,mtrack write(*,*)'track dep # ',ii,' track dep val = ', & dtrack(ii) end do write(*,*)' ' endif c c scoring the crossing info. c if(JTRACK.ne.-1) then call score_cross( & ICPART,IBPART,AMPART,newreg,mreg,xsco,ysco,zsco) endif RETURN * * ======================================================================* * * * Event End DRAWing: * * * *======================================================================* * * ENTRY EEDRAW ( ICODE ) if(idbflg.gt.2) then call dump_common() endif RETURN * *======================================================================* * * * ENergy deposition DRAWing: * * * * Icode = 1x: call from Kaskad * * 10: elastic interaction recoil * * 11: inelastic interaction recoil * * 12: stopping particle * * 13: pseudo-neutron deposition * * 14: escape * * 15: time kill * * Icode = 2x: call from Emfsco * * 20: local energy deposition (i.e. photoelectric) * * 21: below threshold, iarg=1 * * 22: below threshold, iarg=2 * * 23: escape * * 24: time kill * * Icode = 3x: call from Kasneu * * 30: target recoil * * 31: below threshold * * 32: escape * * 33: time kill * * Icode = 4x: call from Kashea * * 40: escape * * 41: time kill * * 42: delta ray stack overflow * * Icode = 5x: call from Kasoph * * 50: optical photon absorption * * 51: escape * * 52: time kill * * * *======================================================================* * * ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO ) * +-------------------------------------------------------------------* c c debug c if (idbflg.lt.0) then IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', & FORM = 'UNFORMATTED' ) WRITE (IODRAW) 0, ICODE, JTRACK, SNGL (ETRACK), & SNGL (WTRACK) WRITE (IODRAW) SNGL (XSCO), SNGL (YSCO), SNGL (ZSCO), & SNGL (RULL) endif endif if(idbflg.gt.1) then write(*,*)'------------- endraw: Ev =',ncase,' -------------' write(*,*)'jtrack = ',jtrack,' mreg = ',mreg write(*,*)'rull = ',rull,' icode = ',icode,' idcurr = ',idcurr write(*,*)'idead(idcurr) = ',idead(idcurr) write(*,*)'x,y,x = ',sngl(XSCO), sngl(YSCO), sngl(ZSCO) WRITE (*,*)'cx,cy,cz = ',SNGL(CXTRCK),SNGL(CYTRCK),SNGL(CZTRCK) endif C C to preserve idead flag correctly C idead (idcurr) = icode C if(idbflg.gt.1) then write(*,*)'ENDRAW: Now idead(idcurr) = ',idead(idcurr) endif C if(jtrack.lt.62) then c************* PARTICLES AND HEAVY IONS WITH JTRACK=-2 IF((JTRACK.GE.1).OR.((JTRACK.LE.-2).AND.(JTRACK.GE.-6)))THEN ICPART = ICHRGE(JTRACK) IBPART = IBARCH(JTRACK) AMPART = AM(JTRACK) c*************HEAVY IONS WITH JTRACK.LT.-6 ELSEIF(JTRACK.LT.-6)THEN ICPART = ICHEAV(-JTRACK) IBPART = IBHEAV(-JTRACK) AMPART = AMNHEA(-JTRACK) ELSE ICPART= ICHRGE(JTRACK) IBPART = IBARCH(JTRACK) AMPART = AM(JTRACK) ENDIF c c if not delta ray below threshold than update the current particle c if( (JTRACK.ne.3).or.(ICODE.ne.22)) then call UpdateCurrPart(mreg,icpart,ibpart,ampart,icode, & xsco,ysco,zsco,JEND) if (jtrack.eq.7) then px(idcurr) = pxf(idcurr) py(idcurr) = pyf(idcurr) pz(idcurr) = pzf(idcurr) endif endif endif if (idbflg.lt.0) then IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)) & , I = 1, MTRACK ), & JBK = 1, NQEMGD ) ENDIF endif if (rull.ne.0) then c c inside the start counter c if( mreg.eq.nregSTC ) then erawSTC = rull equenchedSTC=0 IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) equenchedSTC = dtquen(1,1)*abs_STC END IF call score_STC(mreg,erawSTC,equenchedSTC,xsco, & ysco,zsco,xsco,ysco,zsco) endif c c inside the beam monitor c if( mreg.ge.nregFirstBMN.and.mreg.le.nregLastBMN ) then erawBMN = rull equenchedBMN=0 IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) equenchedBMN = dtquen(1,1)*abs_BMN END IF call score_BMN(mreg,erawBMN,equenchedBMN,xsco, & ysco,zsco,xsco,ysco,zsco) endif c c inside the vertex c if( mreg.le.nregLastVTX .and. mreg.ge.nregFirstVTX ) $ then erawVTX = rull equenchedVTX=0 IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) equenchedVTX = dtquen(1,1)*abs_VTX END IF call score_VTX(mreg,erawVTX,equenchedVTX,xsco, & ysco,zsco,xsco,ysco,zsco) endif c c inside the inner tracker c if( mreg.le.nregLastITR .and. mreg.ge.nregFirstITR ) $ then erawITR = rull equenchedITR=0 IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) equenchedITR = dtquen(1,1)*abs_ITR END IF call score_ITR(mreg,erawITR,equenchedITR,xsco, & ysco,zsco,xsco,ysco,zsco) endif c c inside the second drift chamber c if( mreg.ge.nregFirstMSD .and. mreg.le.nregLastMSD ) then erawMSD = rull equenchedMSD=0 IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) equenchedMSD = dtquen(1,1)*abs_MSD END IF call score_MSD(mreg,erawMSD,equenchedMSD,xsco, & ysco,zsco,xsco,ysco,zsco) endif c c inside the scintillator c if ( mreg.le.nregLastSCN .and. mreg.ge.nregFirstSCN ) then erawSCN = rull equenchedSCN=0 IF ( LQEMGD) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULL, DTQUEN ) equenchedSCN = dtquen(1,1)*abs_SCN END IF call score_SCN(mreg,erawSCN,equenchedSCN,xsco,ysco, & zsco,xsco,ysco,zsco) endif c c inside the calorimeter c if( mreg.le.nregLastCAL .and. mreg.ge.nregFirstCAL )then erawCAL = rull equenchedCAL = 0 IF ( LQEMGD )THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) equenchedCAL = dtquen(1,1)*abs_CAL END IF call score_CAL(mreg,erawCAL,equenchedCAL,xsco,ysco, & zsco,xsco,ysco,zsco) endif endif c c end debug c RETURN * *======================================================================* * * * SOurce particle DRAWing: * * * *======================================================================* * ENTRY SODRAW c if (idbflg.lt.0) then IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', & FORM ='UNFORMATTED' ) END IF WRITE (IODRAW) -NCASE, NPFLKA, NSTMAX, SNGL (TKESUM), & SNGL (WEIPRI) endif if(idbflg.gt.1) then write(*,*)'------------- sodraw: Ev =',ncase,' -------------' endif * +-------------------------------------------------------------------* * | (Radioactive) isotope: it works only for 1 source particle on * | the stack for the time being IF ( ILOFLK (NPFLKA) .GE. 100000 .AND. LRADDC (NPFLKA) ) THEN IARES = MOD ( ILOFLK (NPFLKA), 100000 ) / 100 IZRES = MOD ( ILOFLK (NPFLKA), 10000000 ) / 100000 IISRES = ILOFLK (NPFLKA) / 10000000 IONID = ILOFLK (NPFLKA) if (idbflg.lt.0) then WRITE (IODRAW) ( IONID,SNGL(-TKEFLK(I)), & SNGL (WTFLK(I)), SNGL (XFLK (I)), & SNGL (YFLK (I)), SNGL (ZFLK (I)), & SNGL (TXFLK(I)), SNGL (TYFLK(I)), & SNGL (TZFLK(I)), I = 1, NPFLKA ) endif c * | * +-------------------------------------------------------------------* * | Patch for heavy ions: it works only for 1 source particle on * | the stack for the time being ELSE IF ( ABS (ILOFLK (NPFLKA)) .GE. 10000 ) THEN IONID = ILOFLK (NPFLKA) CALL DCDION ( IONID ) if (idbflg.lt.0) then WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-IONID)), & SNGL (WTFLK(I)), SNGL (XFLK (I)), & SNGL (YFLK (I)), SNGL (ZFLK (I)), & SNGL (TXFLK(I)), SNGL (TYFLK(I)), & SNGL (TZFLK(I)), I = 1, NPFLKA ) endif * | * +-------------------------------------------------------------------* * | Patch for heavy ions: ??? ELSE IF ( ILOFLK (NPFLKA) .LT. -6 ) THEN if (idbflg.lt.0) then WRITE (IODRAW) & ( IONID,SNGL(TKEFLK(I)+AMNHEA(-ILOFLK(NPFLKA))), & SNGL (WTFLK(I)), SNGL (XFLK (I)), & SNGL (YFLK (I)), SNGL (ZFLK (I)), & SNGL (TXFLK(I)), SNGL (TYFLK(I)), & SNGL (TZFLK(I)), I = 1, NPFLKA ) endif * | * +-------------------------------------------------------------------* * | ELSE IONID = ILOFLK (NPFLKA) if (idbflg.lt.0) then WRITE (IODRAW) ( ILOFLK(I), SNGL (TKEFLK(I)+AM(ILOFLK(I))), & SNGL (WTFLK(I)), SNGL (XFLK (I)), & SNGL (YFLK (I)), SNGL (ZFLK (I)), & SNGL (TXFLK(I)), SNGL (TYFLK(I)), & SNGL (TZFLK(I)), I = 1, NPFLKA ) endif END IF c c first call to initialing the kinematic tracking c call UpdateCurrPart(mreg,icpart,ibpart,ampart,icode, & 0.d+00,0.d+00,0.d+00,JSTART) c if(idbflg.gt.0) then write(*,*)'reg = ',numreg(nump),' tprod = ',tempo(nump), & ' jpa = ',jpa(nump) write(*,*)'vert = ',vxi(nump),vyi(nump),vzi(nump), & ' p = ',px(nump),py(nump),pz(nump) endif * | * +-------------------------------------------------------------------* RETURN * *======================================================================* * * * USer dependent DRAWing: * * * * Icode = 99: call from Doiosp, ion splitting secondaries * * Icode = 1xy: call from Kaskad * * 100: elastic interaction secondaries * * 101: inelastic interaction secondaries * * 102: particle decay secondaries * * 103: delta ray generation secondaries * * 104: pair production secondaries * * 105: bremsstrahlung secondaries * * 106: de-excitation in flight secondaries * * 110: radioactive decay products * * Icode = 2xy: call from Emfsco * * 208: bremsstrahlung secondaries * * 210: Moller secondaries * * 212: Bhabha secondaries * * 214: in-flight annihilation secondaries * * 215: annihilation at rest secondaries * * 217: pair production secondaries * * 219: Compton scattering secondaries * * 221: photoelectric secondaries * * 225: Rayleigh scattering secondaries * * 230: inverse Compton secondaries * * 237: mu pair production secondaries * * Icode = 30x: call from Kasneu * * 300: interaction secondaries * * Icode = 40x: call from Kashea * * 400: delta ray generation secondaries * * * * For all interactions secondaries are put on GENSTK common (kp=1,np) * * but for KASHEA delta ray generation where only the secondary elec- * * tron is present and stacked on FLKSTK common for kp=npflka * * * *======================================================================* * ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO ) C C To mark correctly interaction code of current particle C intcode(idcurr) = icode c c debug printing c if (idblfg.lt.0) then IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', & FORM = 'UNFORMATTED' ) END IF endif if(idbflg.gt.0) then call GEOR2N ( mreg, MREGNAM, IERR ) write(*,*)' ' write(*,*)'------------- usdraw: Ev =',ncase,' -------------' write(*,*)' Fluka index= ',ISPUSR(MKBMX2),' idcurr= ',idcurr, & 'jtrack = ',jtrack,' icode = ',icode write(*,*)'regione = ',mregnam, 'x,y,z = ',XSCO,YSCO,ZSCO if(idbflg.gt.2) then write(*,*)'zfrttk = ',zfrttk,' icpart = ',icpart,',m = ', & AMPART write(*,*)'ptrack = ',ptrack,' ekin = ',etrack-AMPART, & ' t = ',atrack WRITE (*,*)'cx,cy,cz = ',SNGL(CXTRCK), SNGL(CYTRCK) & ,SNGL(CZTRCK) endif write(*,*)' USDRAW secondary writing : NP=',NP do ip = 1, NP WRITE(*, *) ip,' jpa sec= ',KPART(ip),' Ek sec= ', & TKI(ip),' p sec= ',plr(ip) end do write(*,*)' USDRAW heavy secondary writing : NPHEAV=',NPHEAV do ip = 1, NPHEAV WRITE(*, *) ip,' jpa sec= ',KHEAVY(ip),' Ek sec= ', & tkheav(ip),' p sec= ',pheavy(ip) end do write(*,*)' ' endif IF (icode.eq.101) THEN !! inelastic interaction IF (MREG.eq.nregSTC) THEN stcfrag = LTRACK if (idbflg.gt.2) then write(*,*) ' mgdraw: STC ',stcfrag do ip = 1, NP if(kpart(ip).lt.-1) then CALL USRDCI(kpart(ip),IONA,IONZ,IONM) write(*,*)'STC: A= ',iona,' Z= ',ionz endif end do do ip = 1, NPHEAV if(kpart(ip).lt.-1) then write(*,*)'STC: A= ',ibheav(ip),' Z= ',icheav(ip) endif end do endif ENDIF IF (MREG.ge.4 .AND. MREG.le.45) THEN bmnfrag = LTRACK if (idbflg.gt.2) then write(*,*) ' mgdraw: BMN ',bmnfrag do ip = 1, NP if(kpart(ip).lt.-1) then CALL USRDCI(kpart(ip),IONA,IONZ,IONM) write(*,*)'BMN: A= ',iona,' Z= ',ionz endif end do do ip = 1, NPHEAV if(kpart(ip).lt.-1) then write(*,*)'BMN: A= ',ibheav(ip),' Z= ',icheav(ip) endif end do endif ENDIF IF (MREG.eq.nregtarg) THEN tgtfrag = LTRACK if (idbflg.gt.2) then write(*,*) ' mgdraw: TGT ',tgtfrag do ip = 1, NP if(kpart(ip).lt.-1) then CALL USRDCI(kpart(ip),IONA,IONZ,IONM) write(*,*)'TGT: A= ',iona,' Z= ',ionz endif end do do ip = 1, NPHEAV if(kpart(ip).lt.-1) then write(*,*)'TGT: A= ',ibheav(ip),' Z= ',icheav(ip) endif end do endif ENDIF IF ( (MREG.eq.nregaria .or. MREG.eq.nregMagAir) ) THEN airfrag = LTRACK if (idbflg.gt.2) then write(*,*) ' mgdraw: AIR ',airfrag do ip = 1, NP if(kpart(ip).lt.-1) then CALL USRDCI(kpart(ip),IONA,IONZ,IONM) write(*,*)'AIR: A= ',iona,' Z= ',ionz endif end do do ip = 1, NPHEAV if(kpart(ip).lt.-1) then write(*,*)'AIR: A= ',ibheav(ip),' Z= ',icheav(ip) endif end do endif ENDIF IF ( MREG.ne.nregaria .AND. MREG.ne.nregMagAir & .AND. MREG.ne.nregSTC & .AND. MREG.ne.nregtarg & .AND. MREG.lt.nregFirstSCN ) THEN elsfrag = LTRACK if (idbflg.gt.2) then write(*,*) ' mgdraw: Else ', elsfrag do ip = 1, NP if(kpart(ip).lt.-1) then CALL USRDCI(kpart(ip),IONA,IONZ,IONM) write(*,*)'AIR: A= ',iona,' Z= ',ionz endif end do do ip = 1, NPHEAV if(kpart(ip).lt.-1) then write(*,*)'AIR: A= ',ibheav(ip),' Z= ',icheav(ip) endif end do endif ENDIF ENDIF * * Do nothing for elastic scattering * +-------------------------------------------------------------------* * | Delta rays, Compton and Moller-Bhabba scattering: IF ( ICODE.EQ.103 .OR. ICODE.EQ.210 .OR. ICODE.EQ. 212 ) THEN ldead = .false. * | * +-------------------------------------------------------------------* * | Bremsstrahlung: ELSE IF ( ICODE .EQ. 208 .OR. ICODE .EQ. 105 ) THEN NPA = NP ldead = .false. * +-------------------------------------------------------------------* * | Compton scattering: ELSE IF ( ICODE .EQ. 219 ) THEN idead (idcurr) = icode ldead = .true. * +-------------------------------------------------------------------* * | PhotoElectric: ELSE IF ( ICODE .EQ. 221 ) THEN idead (idcurr) = icode ldead = .true. * | * +-------------------------------------------------------------------* * | In flight or at rest annihilation: ELSE IF ( ICODE .EQ. 214 .OR. ICODE .EQ. 215 ) THEN idead (idcurr) = icode ldead = .true. * | * +-------------------------------------------------------------------* * | Pair production: ELSE IF ( ICODE.EQ.217 .OR. ICODE.EQ.104 ) THEN idead (idcurr) = icode ldead = .true. * | * +-------------------------------------------------------------------* * | Low energy neutron secondaries: ELSE IF ( ICODE.EQ.300 ) THEN idead (idcurr) = icode ldead =.true. * | * +-------------------------------------------------------------------* * | Decays ELSE IF ( ICODE.EQ.102 .OR. icode.eq.110) THEN idead (idcurr) = icode ldead = .true. * | * +-------------------------------------------------------------------* * | Inelastic interactions: ELSE IF ( ICODE.EQ.101 ) THEN idead (idcurr) = icode ldead = .true. * +-------------------------------------------------------------------* * | In-flight De-excitation: ELSE IF ( ICODE .EQ.106 ) THEN idead (idcurr) = icode ldead = .true. * | * +-------------------------------------------------------------------* * | Everything else: ELSE ldead = .false. END IF * | * +-------------------------------------------------------------------* IF(ICODE.LT.106) THEN numint = numint +1 if (numint.GT.maxint) & write(*,*) 'MGDRAW: numint>maxint',ncase,numint,maxint xint(numint) = xsco yint(numint) = ysco zint(numint) = zsco intpa(numint) = idcurr if(idbflg.gt.1) then write(*,*)'USDRAW: Now numint, idead(numint) = ', & numint,intpa(numint) endif intcode(numint) = intcode(idcurr) ENDIF IF (icode.eq.219 .or. icode.eq.101 .or. icode.eq.102 .or. & icode.eq.110 .or. icode.eq.214 .or. icode.eq.215 .or. & icode.eq.217 .or. icode.eq.221 ) then C call UpdateCurrPart(mreg,0,0,zerzer,icode, & xsco,ysco,zsco,JEND) C endif c if ( ldead ) then idead (idcurr) = icode if(idbflg.gt.1) then write(*,*)'USDRAW: ldead Now idcurr, idead(idcurr) = ', & idcurr,idead(idcurr),icode endif vxf (idcurr) = sngl(xsco) vyf (idcurr) = sngl(ysco) vzf (idcurr) = sngl(zsco) end if c RETURN *=== End of subrutine Mgdraw ==========================================* END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: mgdraw_lib.f C C This is not a standard FLUKA user routine, but it is a service routine for MGDRAW. C It contains the custom subroutines that fill the output arrays for every specific detector and for crossing borders. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c---------------------------------------------------- subroutine dump_common_part() c----------------------------------------------------- c integer ii include 'mgdraw.inc' double precision ptot c write(*,*)' ' write(*,*)'------- dump_common_part ---------' write(*,*)' ' c write(*,*)'DUMP common PART' write(*,*)'no. of particles = ',nump do ii = 1, nump write(*,*)'particle = ', ii write(*,*)'idparent = ', idpa(ii),' igen = ',igen(ii), & ' charge = ',icha(ii),' reg =',numreg(ii) write(*,*)'idead = ',idead(ii),' jpa =',jpa(ii), & ' mass = ',amass(ii),' time = ', tempo(ii) write(*,*)'vert in = ',vxi(ii),vyi(ii),vzi(ii) write(*,*)'vert fin = ',vxf(ii),vyf(ii),vzf(ii) write(*,*)'p ini = ',px(ii),py(ii),pz(ii) ptot = sqrt(px(ii)**2+py(ii)**2+pz(ii)**2) write(*,*)'p tot ini = ',ptot,' ekin ini= ',sqrt(ptot*ptot + & amass(ii)*amass(ii))-amass(ii) write(*,*)'p fin = ',pxf(ii),pyf(ii),pzf(ii) write(*,*)'p tot fin = ',sqrt(pxf(ii)**2+pyf(ii)**2+pzf(ii)**2) end do return end c c c---------------------------------------------------- subroutine dump_common_SCN() c----------------------------------------------------- c integer ii include 'mgdraw.inc' c write(*,*)' ' write(*,*)'------- dump_common_SCN ---------' write(*,*)' ' c write(*,*)'no. of hits in SCN: ',nSCN do ii = 1,nSCN write(*,*)'hit# = ',ii,' part id= ',idSCN(ii) write(*,*)'x,y,z in= ',xinSCN(ii),yinSCN(ii),zinSCN(ii) write(*,*)'x,y,z out= ',xoutSCN(ii),youtSCN(ii), & zoutSCN(ii) write(*,*)'px,py,pz in= ',pxinSCN(ii),pyinSCN(ii), & pzinSCN(ii) write(*,*)'px,py,pz out= ',pxoutSCN(ii),pyoutSCN(ii), & pzoutSCN(ii) write(*,*)'de= ',deSCN(ii),' al= ',alSCN(ii), & ' time= ',timSCN(ii) write(*,*) end do c return end c c c---------------------------------------------------- subroutine dump_common_CROSS() c----------------------------------------------------- c integer ii include 'mgdraw.inc' c write(*,*)' ' write(*,*)'------- dump_common_CROSS ---------' write(*,*)' ' c write(*,*)'no. of crossings: ',nCROSS do ii = 1,nCROSS write(*,*)'crossing# = ',ii,' part id = ',idCROSS(ii) write(*,*)'x,y,z = ',xCROSS(ii),yCROSS(ii),zCROSS(ii) write(*,*)'px,py,pz in = ',pxCROSS(ii),pyCROSS(ii), & pzCROSS(ii) write(*,*)'reg = ',nregCROSS(ii),' regold = ',nregoldCROSS(ii), & ' m = ',amaCROSS(ii),' time = ',tCROSS(ii), & ' cha= ',chCROSS(ii) write(*,*) end do c return end c c c---------------------------------------------------- subroutine dump_common() c----------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' include 'mgdraw.inc' c if(idbflg .gt. 0) then c call dump_common_part() c call dump_common_SCN() c call dump_common_CROSS() c endif c return end c c c------------------------------------------------------------------------- subroutine score_STC(mreg,erawSTC,equenchedSTC,xcordin, & ycordin,zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawSTC, equenchedSTC include 'mgdraw.inc' integer ii, nSTC_now c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_STC-----------------' write(*,*)'idcurr= ',idcurr,' equenchedSTC= ', & equenchedSTC write(*,*)'dtrack= ',(dtrack(ii),ii=1,mtrack) endif c c new hit in the start counter? c nSTC_now = 0 do ii=1,nSTC if(idSTC(ii).eq.idcurr) nSTC_now = ii end do c c if nSTC_now = 0 then ... new hit... c if(nSTC_now.eq.0) then if(nSTC.eq.maxSTC) then write(*,*)'Error: Score_STC:' write(*,*)'Maximum number of STC scoring exceeded: ev= ', & ncase return endif nSTC = nSTC + 1 nSTC_now = nSTC idSTC(nSTC) = idcurr xinSTC(nSTC_now) = sngl(xcordin) yinSTC(nSTC_now) = sngl(ycordin) zinSTC(nSTC_now) = sngl(zcordin) pxinSTC(nSTC_now) = sngl(ptrack*cxtrck) pyinSTC(nSTC_now) = sngl(ptrack*cytrck) pzinSTC(nSTC_now) = sngl(ptrack*cztrck) timSTC(nSTC_now) = sngl(atrack) deSTC(nSTC_now) = 0. endif c xoutSTC(nSTC_now) = sngl(xcordout) youtSTC(nSTC_now) = sngl(ycordout) zoutSTC(nSTC_now) = sngl(zcordout) pxoutSTC(nSTC_now) = sngl(ptrack*cxtrck) pyoutSTC(nSTC_now) = sngl(ptrack*cytrck) pzoutSTC(nSTC_now) = sngl(ptrack*cztrck) deSTC(nSTC_now) = deSTC(nSTC_now) + sngl(erawSTC) alSTC(nSTC_now)=alSTC(nSTC_now)+ sngl(equenchedSTC) c if(idbflg.gt.1) then write(*,*)'nSTC= ',nSTC_now, & ' deSTC(nSTC)= ',deSTC(nSTC_now) write(*,*)' ' endif c return end c c c------------------------------------------------------------------------- subroutine score_BMN(mreg,erawBMN,equenchedBMN, & xcordin,ycordin,zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawBMN, equenchedBMN include 'mgdraw.inc' include 'parameters.inc' integer ii, nBMN_now integer iview, ilay, icell c iview = ireg2viewBMN(mreg) ilay = ireg2layBMN(mreg) icell = ireg2cellBMN(mreg) c if ((ilay.ge.nlayBMN.or.ilay.lt.0).or. & (icell.ge.ncellBMN.or.icell.lt.0).or. & (iview.ne.1.and.iview.ne.-1)) then write(*,*) ' WARNING!!!! ilay, iview, icell, mreg ' write(*,*) ilay, iview, icell, mreg endif c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_BMN-----------------' write(*,*)'idcurr= ',idcurr,' equenchedBMN= ', & equenchedBMN,' view= ', iview,' lay= ', ilay, & ' cell= ',icell,' mreg= ',mreg write(*,*)'xyz= ',xcordin,ycordin,zcordin write(*,*)'dtrack= ',(dtrack(ii),ii=1,mtrack) endif c c new hit in the beam monitor? c nBMN_now = 0 do ii=1,nBMN if( (idBMN(ii).eq.idcurr) .and. (iviewBMN(ii).eq.iview) .and. & (ilayBMN(ii).eq.ilay) .and. (icellBMN(ii).eq.icell) )then nBMN_now = ii endif end do c c if nBMN_now = 0 then ... new hit... c if(nBMN_now.eq.0) then if(nBMN.eq.maxBMN) then write(*,*)'Error: Score_BMN:' write(*,*)'Maximum number of BMN scoring exceeded: ev= ', & ncase return endif nBMN = nBMN + 1 nBMN_now = nBMN idBMN(nBMN) = idcurr xinBMN(nBMN_now) = sngl(xcordin) yinBMN(nBMN_now) = sngl(ycordin) zinBMN(nBMN_now) = sngl(zcordin) pxinBMN(nBMN_now) = sngl(ptrack*cxtrck) pyinBMN(nBMN_now) = sngl(ptrack*cytrck) pzinBMN(nBMN_now) = sngl(ptrack*cztrck) timBMN(nBMN_now) = sngl(atrack) iviewBMN(nBMN_now) = iview ilayBMN(nBMN_now) = ilay icellBMN(nBMN_now) = icell deBMN(nBMN_now) = 0. endif c xoutBMN(nBMN_now) = sngl(xcordout) youtBMN(nBMN_now) = sngl(ycordout) zoutBMN(nBMN_now) = sngl(zcordout) pxoutBMN(nBMN_now) = sngl(ptrack*cxtrck) pyoutBMN(nBMN_now) = sngl(ptrack*cytrck) pzoutBMN(nBMN_now) = sngl(ptrack*cztrck) deBMN(nBMN_now) = deBMN(nBMN_now) + sngl(erawBMN) alBMN(nBMN_now)=alBMN(nBMN_now)+ sngl(equenchedBMN) c if(idbflg.gt.1) then write(*,*)'nBMN= ',nBMN_now, ' deBMN(nBMN)= ',deBMN(nBMN_now) write(*,*)' ' endif c return end c c------------------------------------------------------------------------- subroutine score_VTX(mreg, erawVTX, equenchedVTX,xcordin, & ycordin,zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawVTX, equenchedVTX include 'mgdraw.inc' include 'parameters.inc' integer mreg, ii integer ilay c ilay = mreg-nregFirstVTX if (ilay.ge.nlayVTX.or.ilay.lt.0) THEN write(*,*) ' WARNING!!!! ilayVTX= ',ilay, & ' zcordin= ',zcordin ENDIF c if(ilay.lt.0.or.ilay.ge.nlayVTX)then write(*,*)'WARNING! max vertex lay no. exceeded: ilay= ', & ilay,' mreg= ',mreg endif c nVTX_now = 0 c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_VTX-----------------' write(*,*)'test ','idcurr= ',idcurr,' lay= ', & ilay, ' mreg= ',mreg,' nVTX= ',nVTX write(*,*)'xyz= ',xcordin,ycordin,zcordin endif c c new hit in the vertex? c if (idbflg.gt.1) then write(*,*)'Check of previous hits nVTX_now= ',nVTX_now, & ' nVTX= ',nVTX endif do ii=1,nVTX if (idbflg.gt.1) then write(*,*)'ii= ',ii write(*,*)'idVTX(ii)= ',idVTX(ii),' idcurr= ',idcurr write(*,*)'ilay= ',ilay endif if( ( idVTX(ii).eq.idcurr .and. & ilayVTX(ii).eq.ilay) .OR. & ( sngl(xcordin).eq.xinVTX(ii) .and. & sngl(ycordin).eq.yinVTX(ii) .and. & sngl(zcordin).eq.zinVTX(ii) ) & ) then nVTX_now = ii if (idbflg.gt.1) then write(*,*)'Previous hit found: nVTX_now= ',nVTX_now, & ' ii= ',ii endif endif end do c c if nVTX_now = 0 then ... new hit... c if( nVTX_now.eq.0) then if(nVTX.eq.maxVTX) then write(*,*)'Error: Score_VTX:' write(*,*)'Maximum hit number exceeded : ev= ',ncase return endif nVTX = nVTX + 1 nVTX_now = nVTX idVTX(nVTX) = idcurr ilayVTX(nVTX_now) = ilay xinVTX(nVTX_now) = sngl(xcordin) yinVTX(nVTX_now) = sngl(ycordin) zinVTX(nVTX_now) = sngl(zcordin) pxinVTX(nVTX_now) = sngl(ptrack*cxtrck) pyinVTX(nVTX_now) = sngl(ptrack*cytrck) pzinVTX(nVTX_now) = sngl(ptrack*cztrck) timVTX(nVTX_now) = sngl(atrack) deVTX(nVTX_now) = 0. endif c xoutVTX(nVTX_now) = sngl(xcordout) youtVTX(nVTX_now) = sngl(ycordout) zoutVTX(nVTX_now) = sngl(zcordout) pxoutVTX(nVTX_now) = sngl(ptrack*cxtrck) pyoutVTX(nVTX_now) = sngl(ptrack*cytrck) pzoutVTX(nVTX_now) = sngl(ptrack*cztrck) deVTX(nVTX_now) = deVTX(nVTX_now) + sngl(erawVTX) alVTX(nVTX_now) = alVTX(nVTX_now) + sngl(equenchedVTX) c if(idbflg.gt.1) then do ii=1,nVTX write(*,*)'ii= ',ii write(*,*)'idVTX(ii)= ',idVTX(ii),' idcurr= ',idcurr write(*,*)'ilay= ',ilayVTX(ii) enddo write(*,*)'nVTX_now= ',nVTX_now,' nVTX= ',nVTX, & ' deVTX(nVTX)= ',deVTX(nVTX_now) write(*,*)' ' endif c return end c c c------------------------------------------------------------------------- subroutine score_ITR(mreg, erawITR, equenchedITR,xcordin, & ycordin,zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawITR, equenchedITR double precision xmin, xmax, ymin, ymax include 'mgdraw.inc' include 'parameters.inc' integer mreg, ii integer isens isens = ireg2sensITR(mreg) if (isens.ge.nsensITR.or.sens.lt.0) THEN write(*,*) ' WARNING!!!! isensITR= ',isens, & ' xcordin= ',xcordin, ' ycordin= ',ycordin, & ' zcordin= ',zcordin,' mreg= ',mreg ENDIF nITR_now = 0 c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_ITR-----------------' write(*,*)'test ','idcurr= ',idcurr, & ' sens= ', isens, c & ' lay= ', ilay, ' plume=',iplume, ' mimosa= ',imimo, & ' mreg =',mreg, ' nITR= ',nITR write(*,*)'xyz= ',xcordin,ycordin,zcordin endif c c new hit in the inner tracker? c if (idbflg.gt.1) then write(*,*)'Check of previous hits nITR_now= ',nITR_now, & ' nITR= ',nITR endif do ii=1,nITR if (idbflg.gt.1) then write(*,*)'ii= ',ii write(*,*)'idITR(ii)= ',idITR(ii),' idcurr= ',idcurr write(*,*)'isens= ',isens c write(*,*)'iplume= ',iplume,'ilay= ',ilay,'imimo= ',imimo endif if( ( idITR(ii).eq.idcurr .and. & isensITR(ii).eq.isens ) .OR. & ( sngl(xcordin).eq.xinITR(ii) .and. & sngl(ycordin).eq.yinITR(ii) .and. & sngl(zcordin).eq.zinITR(ii) ) & ) then nITR_now = ii if (idbflg.gt.1) then write(*,*)'Previous hit found: nITR_now= ',nITR_now, & ' ii= ',ii endif endif end do c c if nITR_now = 0 then ... new hit... c if( nITR_now.eq.0) then if(nITR.eq.maxITR) then write(*,*)'Error: Score_ITR:' write(*,*)'Maximum hit number exceeded : ev= ',ncase return endif nITR = nITR + 1 nITR_now = nITR idITR(nITR) = idcurr isensITR(nITR_now) = isens xinITR(nITR_now) = sngl(xcordin) yinITR(nITR_now) = sngl(ycordin) zinITR(nITR_now) = sngl(zcordin) pxinITR(nITR_now) = sngl(ptrack*cxtrck) pyinITR(nITR_now) = sngl(ptrack*cytrck) pzinITR(nITR_now) = sngl(ptrack*cztrck) timITR(nITR_now) = sngl(atrack) deITR(nITR_now) = 0. endif c xoutITR(nITR_now) = sngl(xcordout) youtITR(nITR_now) = sngl(ycordout) zoutITR(nITR_now) = sngl(zcordout) pxoutITR(nITR_now) = sngl(ptrack*cxtrck) pyoutITR(nITR_now) = sngl(ptrack*cytrck) pzoutITR(nITR_now) = sngl(ptrack*cztrck) deITR(nITR_now) = deITR(nITR_now) + sngl(erawITR) alITR(nITR_now) = alITR(nITR_now) + sngl(equenchedITR) c if(idbflg.gt.1) then do ii=1,nITR write(*,*)'ii= ',ii write(*,*)'idITR(ii)= ',idITR(ii),' idcurr= ',idcurr c write(*,*)'ilay= ',ilayITR(ii) write(*,*)'isens= ',isensITR(ii) enddo write(*,*)'nITR_now= ',nITR_now,' nITR= ',nITR, & ' deITR(nITR)= ',deITR(nITR_now) write(*,*)' ' endif c return end c c c------------------------------------------------------------------------- subroutine score_MSD(mreg,erawMSD,equenchedMSD,xcordin, & ycordin,zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawMSD, equenchedMSD include 'mgdraw.inc' include 'parameters.inc' integer ii, nMSD_now integer ilay c ilay = ireg2layMSD(mreg) c if (ilay.ge.nlayMSD.or.ilay.lt.0) then write(*,*) ' WARNING!!!! ilay ', ilay endif c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_MSD-----------------' write(*,*)'idcurr= ',idcurr,' equenchedMSD= ', & equenchedMSD,' lay= ', ilay, & ' cell= ',icell,' mreg= ',mreg write(*,*)'xyz= ',xcordin,ycordin,zcordin write(*,*)'dtrack= ',(dtrack(ii),ii=1,mtrack) endif c c new hit in the second drift chamber? c nMSD_now = 0 do ii=1,nMSD if((idMSD(ii).eq.idcurr).and. & (ilayMSD(ii).eq.ilay))then nMSD_now = ii endif end do c c if nMSD_now = 0 then ... new hit... c if(nMSD_now.eq.0) then if(nMSD.eq.maxMSD) then write(*,*)'Error: Score_MSD:' write(*,*)'Maximum number of MSD scoring exceeded: ev= ', & ncase return endif nMSD = nMSD + 1 nMSD_now = nMSD idMSD(nMSD) = idcurr xinMSD(nMSD_now) = sngl(xcordin) yinMSD(nMSD_now) = sngl(ycordin) zinMSD(nMSD_now) = sngl(zcordin) pxinMSD(nMSD_now) = sngl(ptrack*cxtrck) pyinMSD(nMSD_now) = sngl(ptrack*cytrck) pzinMSD(nMSD_now) = sngl(ptrack*cztrck) timMSD(nMSD_now) = sngl(atrack) ilayMSD(nMSD_now) = ilay deMSD(nMSD_now) = 0. endif c xoutMSD(nMSD_now) = sngl(xcordout) youtMSD(nMSD_now) = sngl(ycordout) zoutMSD(nMSD_now) = sngl(zcordout) pxoutMSD(nMSD_now) = sngl(ptrack*cxtrck) pyoutMSD(nMSD_now) = sngl(ptrack*cytrck) pzoutMSD(nMSD_now) = sngl(ptrack*cztrck) deMSD(nMSD_now) = deMSD(nMSD_now) + sngl(erawMSD) alMSD(nMSD_now)=alMSD(nMSD_now)+ sngl(equenchedMSD) c if(idbflg.gt.1) then write(*,*)'nMSD= ',nMSD_now, & ' deMSD(nMSD)= ',deMSD(nMSD_now) write(*,*)' ' endif c return end c c c------------------------------------------------------------------------- subroutine score_SCN(mreg,erawSCN,equenchedSCN,xcordin, & ycordin,zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawSCN, equenchedSCN include 'mgdraw.inc' include 'parameters.inc' integer ii, nSCN_now integer stripSCN, viewSCN c stripSCN = ireg2stripSCN(mreg) viewSCN = ireg2viewSCN(mreg) if (stripSCN.ge.nstripSCN .or. colSCN.lt.0) THEN write(*,*) ' WARNING!!!! stripSCN = ', stripSCN, & ' xcordin= ',xcordin, ' ycordin= ',ycordin, & ' zcordin= ',zcordin ENDIF if (viewSCN.gt.1 .or. viewSCN.lt.-1) THEN write(*,*) ' WARNING!!!! viewSCN = ', viewSCN, & ' xcordin= ',xcordin, ' ycordin= ',ycordin, & ' zcordin= ',zcordin ENDIF c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_SCN-----------------' write(*,*)'idcurr= ',idcurr,' equenchedSCN= ', & equenchedSCN write(*,*)'dtrack= ',(dtrack(ii),ii=1,mtrack) write(*,*)'mreg= ',mreg,' strip= ',stripSCN, & ' view= ',viewSCN write(*,*)'xyz= ',xcordin,ycordin,zcordin endif c c new hit in the scintillator? c nSCN_now = 0 do ii=1,nSCN if((idSCN(ii).eq.idcurr.and. & iviewSCN(ii).eq.viewSCN.and. & istripSCN(ii).eq.stripSCN ) .OR. & ( sngl(xcordin).eq.xinSCN(ii) .and. & sngl(ycordin).eq.yinSCN(ii) .and. & sngl(zcordin).eq.zinSCN(ii) ))then nSCN_now = ii endif end do c c if nSCN_now = 0 then ... new hit... c if(nSCN_now.eq.0) then if(nSCN.eq.maxSCN ) then write(*,*)'Error: Score_SCN:' write(*,*)'Maximum number of SCN scoring exceeded : ev= ', & ncase return endif nSCN = nSCN + 1 nSCN_now = nSCN idSCN(nSCN) = idcurr istripSCN(nSCN_now) = stripSCN iviewSCN(nSCN_now) = viewSCN xinSCN(nSCN_now) = sngl(xcordin) yinSCN(nSCN_now) = sngl(ycordin) zinSCN(nSCN_now) = sngl(zcordin) pxinSCN(nSCN_now) = sngl(ptrack*cxtrck) pyinSCN(nSCN_now) = sngl(ptrack*cytrck) pzinSCN(nSCN_now) = sngl(ptrack*cztrck) timSCN(nSCN_now) = sngl(atrack) deSCN(nSCN_now) = 0. endif c xoutSCN(nSCN_now) = sngl(xcordout) youtSCN(nSCN_now) = sngl(ycordout) zoutSCN(nSCN_now) = sngl(zcordout) pxoutSCN(nSCN_now) = sngl(ptrack*cxtrck) pyoutSCN(nSCN_now) = sngl(ptrack*cytrck) pzoutSCN(nSCN_now) = sngl(ptrack*cztrck) deSCN(nSCN_now) = deSCN(nSCN_now) + sngl(erawSCN) alSCN(nSCN_now) = alSCN(nSCN_now) + sngl(equenchedSCN) c if(idbflg.gt.1) then write(*,*)'nSCN= ',nSCN_now, & ' deSCN(nSCN)= ',deSCN(nSCN_now) write(*,*)'istripSCN= ',istripSCN(nSCN_now), & 'iviewSCN= ',iviewSCN(nSCN_now) write(*,*)' ' endif c return end c c c------------------------------------------------------------------------- subroutine score_CAL(mreg, erawCAL, equenchedCAL,xcordin,ycordin, & zcordin,xcordout,ycordout,zcordout) c-------------------------------------------------------------------------- INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' double precision erawCAL, equenchedCAL include 'mgdraw.inc' include 'parameters.inc' integer ii, mreg, nCAL_now integer cryCAL c cryCAL = ireg2cryCAL(mreg) if (cryCAL.ge.ncryCAL .or. cryCAL.lt.0) THEN write(*,*) ' WARNING!!!! cryCAL = ', cryCAL, & ' xcordin= ',xcordin, ' ycordin= ',ycordin ENDIF c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'---------------Score_CAL-----------------' write(*,*)'idcurr= ',idcurr,' equenchedlyso= ', & equenchedCAL write(*,*)'dtrack= ',(dtrack(ii),ii=1,mtrack) endif c nCAL_now = 0 c c new hit in the calo? c do ii=1,nCAL if( ((idCAL(ii).eq.idcurr) .and. (icryCAL(ii).eq.cryCAL)) .or. & ( sngl(xcordin).eq.xinCAL(ii) .and. & sngl(ycordin).eq.yinCAL(ii) .and. & sngl(zcordin).eq.zinCAL(ii) )) then nCAL_now = ii endif end do c c if nCAL_now = 0 then ... new hit... c if(nCAL_now.eq.0) then if(nCAL.eq.maxCAL) then write(*,*)'Error: Score_CAL:' write(*,*)'Maximum number of CAL scorings exceeded : ev= ', & ncase return endif nCAL = nCAL + 1 nCAL_now = nCAL idCAL(nCAL_now) = idcurr icryCAL(nCAL_now) = cryCAL xinCAL(nCAL_now) = sngl(xcordin) yinCAL(nCAL_now) = sngl(ycordin) zinCAL(nCAL_now) = sngl(zcordin) pxinCAL(nCAL_now) = sngl(ptrack*cxtrck) pyinCAL(nCAL_now) = sngl(ptrack*cytrck) pzinCAL(nCAL_now) = sngl(ptrack*cztrck) timCAL(nCAL_now) = sngl(atrack) deCAL(nCAL_now) = 0. endif c xoutCAL(nCAL_now) = sngl(xcordout) youtCAL(nCAL_now) = sngl(ycordout) zoutCAL(nCAL_now) = sngl(zcordout) pxoutCAL(nCAL_now) = sngl(ptrack*cxtrck) pyoutCAL(nCAL_now) = sngl(ptrack*cytrck) pzoutCAL(nCAL_now) = sngl(ptrack*cztrck) deCAL(nCAL_now) = deCAL(nCAL_now) + sngl(erawCAL) alCAL(nCAL_now) = alCAL(nCAL_now) + sngl(equenchedCAL) c if(idbflg.gt.0) then write(*,*)' ' write(*,*)'---------------Score_CAL-----------------' write(*,*)'idcurr= ',idcurr,' equenchedlyso= ', & equenchedCAL write(*,*)'dtrack= ',(dtrack(ii),ii=1,mtrack) endif if(idbflg.gt.1) then write(*,*)'nCAL= ',nCAL_now, & ' deCAL(nCAL)= ',deCAL(nCAL_now) write(*,*)'icryCAL= ',icryCAL(nCAL_now) write(*,*)' ' endif c return end c c c------------------------------------------------------------------------- SUBROUTINE score_CROSS( & icharge,numbar,ampart,newreg,mreg,xsco,ysco,zsco) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' include 'mgdraw.inc' integer icharge, numbar, newreg, mreg double precision ampart, xsco, ysco, zsco c if(nCROSS.ge.maxCROSS) then write(*,*)"SCORE_CROSS: max number of CROSSing exceeded" return endif nCROSS = nCROSS + 1 idCROSS(nCROSS) = idcurr nregCROSS(nCROSS) = newreg nregoldCROSS(nCROSS) = mreg xCROSS(nCROSS) = sngl(xsco) yCROSS(nCROSS) = sngl(ysco) zCROSS(nCROSS) = sngl(zsco) pxCROSS(nCROSS) = sngl(ptrack*cxtrck) pyCROSS(nCROSS) = sngl(ptrack*cytrck) pzCROSS(nCROSS) = sngl(ptrack*cztrck) tCROSS(nCROSS) = sngl(atrack) chCROSS(nCROSS) = icharge amaCROSS(nCROSS) = sngl(ampart) c if(idbflg.gt.1) then write(*,*)' ' write(*,*)'--------------- Score_CROSS -----------------' write(*,*)'idcurr = ',idcurr,' nCROSS= ',nCROSS write(*,*)'reg= ',newreg,' mreg= ',mreg, & ' pxy,z= ',pxCROSS(nCROSS), & pyCROSS(nCROSS),pzCROSS(nCROSS),' mass= ',ampart write(*,*)'x,y,zCROSS = ',xsco,ysco,zsco,' t= ',atrack, & ' cha= ',icharge endif c return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: usreou.f C C This subroutine is called at the end of each event, i.e. when the primary and all its descendant particles have been C transported. In the FOOT case, this subroutine writes the output arrays on the output ASCII file at the end of the event. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE USREOU INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1991-2005 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * USeR Event OUtput: this routine is called at the end of each * * event * * * * Created on 01 january 1991 by Alfredo Ferrari & Paola Sala * * Infn - Milan * * * * Last change on 09-apr-99 by Alfredo Ferrari * * * * * *----------------------------------------------------------------------* * include '(CASLIM)' include '../mgdraw.inc' integer ii logical trigger double precision Eqcalo, Eqscint c c C idbflg = 1 trigger = .false. tarfrag = stcfrag + 10*bmnfrag + 100*tgtfrag + 1000*airfrag & + 10000*elsfrag if (stcfrag.eq.1) trig_stc = .true. if (bmnfrag.eq.1) trig_bmn = .true. if (tgtfrag.eq.1) trig_tgt = .true. if (airfrag.eq.1) trig_air = .true. if (elsfrag.eq.1) trig_els = .true. if(idbflg.gt.0) then write(*,*) '---------------------------' write(*,*)'Usreou. Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag write(*,*) endif Eqcalo = 0 Eqscint = 0 c do ii =1,nCAL Eqcalo = Eqcalo + deCAL(ii) end do c do ii =1,nSCN Eqscint = Eqscint + deSCN(ii) end do C C fragtrig = 0 ----- no trigger C fragtrig = 1 ----- only STC with inelastic interaction C fragtrig = 2 ----- only BMN with inelastic interaction C fragtrig = 3 ----- only TGT with inelastic interaction C fragtrig = 4 ----- only AIR with inelastic interaction C fragtrig = 5 ----- inelastic interaction in STC + X C fragtrig = 6 ----- inelastic interaction in TGT + X C fragtrig = 7 ----- inelastic interaction in AIR + X C fragtrig = 8 ----- inelastic interaction in ELS + X C fragtrig = 9 ----- inelastic interaction in BMN + X C fragtrig = 10 ---- inelastic interaction in STC - BMN +X C fragtrig = 11 ---- inelastic interaction in BMN - STC +X C fragtrig = 34 ---- inelastic interaction in TGT + AIR - STC -BMN -ELS C fragtrig = 35 ---- inelastic interaction in TGT + ELS - STC -BMN -AIR C fragtrig = -1 ---- no inelastic interaction in STC C fragtrig = -2 ---- no inelastic interaction in BMN C fragtrig = -3 ---- no inelastic interaction in TGT C fragtrig = -4 ---- no inelastic interaction in AIR C fragtrig = -123 -- no inelastic interaction in STC, BMN and TGT C fragtrig = -1234 - no inelastic interaction in STC, BMN, TGT and AIR C if (fragtrig.eq.0 ) then trigger = .true. else if (fragtrig.eq.1) then if (trig_stc .and. .not.trig_bmn .and. .not.trig_tgt .and. & .not.trig_air .and. .not.trig_els) then trigger = .true. if(idbflg.gt.0) then write(*,*)'STC Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.2) then if (trig_bmn .and. .not.trig_stc .and. .not.trig_tgt .and. & .not.trig_air .and. .not.trig_els) then trigger = .true. if(idbflg.gt.0) then write(*,*)'BMN Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.3) then if (trig_tgt .and. .not.trig_stc .and. .not.trig_bmn .and. & .not.trig_air .and. .not.trig_els) then trigger = .true. if(idbflg.gt.0) then write(*,*)'TGT Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.4) then if (trig_air .and. .not.trig_stc .and. .not.trig_bmn .and. & .not.trig_tgt .and. .not.trig_els) then trigger = .true. if(idbflg.gt.0) then write(*,*)'AIR Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.5) then if (trig_stc) then trigger = .true. if(idbflg.gt.0) then write(*,*)'STC+X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.6) then if (trig_tgt) then trigger = .true. if(idbflg.gt.0) then write(*,*)'TGT+X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.7) then if (trig_air) then trigger = .true. if(idbflg.gt.0) then write(*,*)'AIR+X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.8) then if (trig_els) then trigger = .true. if(idbflg.gt.0) then write(*,*)'ELS+X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.9) then if (trig_bmn) then trigger = .true. if(idbflg.gt.0) then write(*,*)'BMN+X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.10) then if (trig_stc .and. .not.trig_bmn) then trigger = .true. if(idbflg.gt.0) then write(*,*)'SCT -BMN +X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.11) then if (trig_bmn .and. .not.trig_stc) then trigger = .true. if(idbflg.gt.0) then write(*,*)'BMN -STC +X Ev= ',numev, & ' trig_stc: ', trig_stc, & ' trig_bmn: ', trig_bmn, & ' trig_tgt: ', trig_tgt, & ' trig_air: ', trig_air, & ' trig_els: ', trig_els, & ' trigger= ', Trigger, & ' fragtrig= ',fragtrig, & ' tarfrag= ',tarfrag endif endif else if (fragtrig.eq.-1) then if (.not.trig_stc) then trigger = .true. endif else if (fragtrig.eq.-2) then if (.not.trig_bmn) then trigger = .true. endif else if (fragtrig.eq.-3) then if (.not.trig_tgt) then trigger = .true. endif else if (fragtrig.eq.-4) then if (.not.trig_air) then trigger = .true. endif else if (fragtrig.eq.-5) then if (.not.trig_els) then trigger = .true. endif else if (fragtrig.eq.-123) then if (.not.trig_stc .and. .not.trig_bmn & .and. .not.trig_tgt) then trigger = .true. endif else if (fragtrig.eq.-1234) then if (.not.trig_stc .and. .not.trig_bmn & .and. .not.trig_tgt .and. .not.trig_air) then trigger = .true. endif else if (fragtrig.eq.35) then if (trig_tgt .and. trig_els .and. .not.trig_stc & .and. .not.trig_bmn .and. .not.trig_air) then trigger = .true. endif else if (fragtrig.eq.34) then if (trig_tgt .and. trig_air .and. .not.trig_stc & .and. .not.trig_bmn .and. .not.trig_els) then trigger = .true. endif endif c C idbflg = 0 if(idbflg.gt.0) then write(*,*)'ev= ',numev,' Equench= ',Eqcalo, & ' Eqscint= ',Eqscint, & ' trigger= ', Trigger, & ' Eqcalo= ',Eqcalo,' Eqscint= ',Eqscint, & ' Ethrdep= ',Ethrdep,' tarfrag= ',tarfrag endif c c if (idbflg.eq.10) trigger = .false. if(trigger) then c write(outunit,*) ncase,nump,nSTC,nBMN,nVTX,nITR,nMSD,nSCN,nCAL, & nCROSS c c scrivo la banca delle particelle c do ii = 1,nump write(outunit,*)idpa(ii), igen(ii), icha(ii), & numreg(ii), iba(ii), idead(ii), jpa(ii), vxi(ii), & vyi(ii), vzi(ii), vxf(ii), vyf(ii), vzf(ii), px(ii), & py(ii),pz(ii),pxf(ii),pyf(ii),pzf(ii),amass(ii), & tempo(ii),tof(ii),trlen(ii) end do c c scrivo lo start counter c do ii = 1,nSTC write(outunit,*) idSTC(ii), & xinSTC(ii), yinSTC(ii), zinSTC(ii), & xoutSTC(ii), youtSTC(ii), zoutSTC(ii), & pxinSTC(ii), pyinSTC(ii), pzinSTC(ii), & pxoutSTC(ii), pyoutSTC(ii), pzoutSTC(ii), & deSTC(ii), alSTC(ii), timSTC(ii) end do c c scrivo il beam monitor c do ii = 1,nBMN write(outunit,*) idBMN(ii), iviewBMN(ii), ilayBMN(ii), & icellBMN(ii), xinBMN(ii), yinBMN(ii), zinBMN(ii), & xoutBMN(ii), youtBMN(ii), zoutBMN(ii), & pxinBMN(ii), pyinBMN(ii), pzinBMN(ii), & pxoutBMN(ii), pyoutBMN(ii), pzoutBMN(ii), & deBMN(ii), alBMN(ii), timBMN(ii) end do c c scrivo il vertice c do ii = 1,nVTX write(outunit,*) idVTX(ii), ilayVTX(ii), & xinVTX(ii), yinVTX(ii), zinVTX(ii), & xoutVTX(ii), youtVTX(ii), zoutVTX(ii), & pxinVTX(ii), pyinVTX(ii), pzinVTX(ii), & pxoutVTX(ii), pyoutVTX(ii), pzoutVTX(ii), & deVTX(ii), alVTX(ii), timVTX(ii) end do c c scrivo l'inner tracker c do ii = 1,nITR write(outunit,*) idITR(ii), isensITR(ii), & xinITR(ii), yinITR(ii), zinITR(ii), & xoutITR(ii), youtITR(ii), zoutITR(ii), & pxinITR(ii), pyinITR(ii), pzinITR(ii), & pxoutITR(ii), pyoutITR(ii), pzoutITR(ii), & deITR(ii), alITR(ii), timITR(ii) end do c c scrivo le microstrip c do ii = 1,nMSD write(outunit,*) idMSD(ii), & ilayMSD(ii), & xinMSD(ii), yinMSD(ii), zinMSD(ii), & xoutMSD(ii), youtMSD(ii), zoutMSD(ii), & pxinMSD(ii), pyinMSD(ii), pzinMSD(ii), & pxoutMSD(ii), pyoutMSD(ii), pzoutMSD(ii), & deMSD(ii), alMSD(ii), timMSD(ii) end do c c scrivo lo scintillatore c do ii = 1,nSCN write(outunit,*) idSCN(ii), & istripSCN(ii), iviewSCN(ii), & xinSCN(ii), yinSCN(ii), zinSCN(ii), & xoutSCN(ii), youtSCN(ii), zoutSCN(ii), & pxinSCN(ii), pyinSCN(ii), pzinSCN(ii), & pxoutSCN(ii), pyoutSCN(ii), pzoutSCN(ii), & deSCN(ii), alSCN(ii), timSCN(ii) end do c c scrivo il calorimetro c do ii = 1,nCAL write(outunit,*) idCAL(ii), icryCAL(ii), & xinCAL(ii), yinCAL(ii), zinCAL(ii), & xoutCAL(ii), youtCAL(ii), zoutCAL(ii), & pxinCAL(ii), pyinCAL(ii), pzinCAL(ii), & pxoutCAL(ii), pyoutCAL(ii), pzoutCAL(ii), & deCAL(ii), alCAL(ii), timCAL(ii) end do c c scrivo i crossings c do ii = 1,nCROSS write(outunit,*) idCROSS(ii), nregCROSS(ii), & nregoldCROSS(ii), xCROSS(ii), yCROSS(ii), zCROSS(ii), & pxCROSS(ii), pyCROSS(ii), pzCROSS(ii), & amaCROSS(ii), chCROSS(ii), tCROSS(ii) end do c c c endif c if(idbflg .gt. 0) then call dump_common() write(*,*)'***************************************************' write(*,*)'***************** EVENT END ***********************' write(*,*)'***************************************************' write(*,*)'' endif c c c RETURN *=== End of subroutine Usreou =========================================* END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: usrout.f C C It is called at the end of the run, and is activated by the USROCALL card in the input file. C It can be used to print customized output. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE USROUT ( WHAT, SDUM ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1991-2005 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * USeR OUTput: this routine is called every time the USROCALL card * * is found in the input stream * * * * * * Created on 01 january 1991 by Alfredo Ferrari & Paola Sala * * Infn - Milan * * * * Last change on 20-mar-05 by Alfredo Ferrari * * * * * *----------------------------------------------------------------------* * DIMENSION WHAT (6) CHARACTER SDUM*8 c ITSAV = ITEXPI NLUNI1 = 63 CLOSE(NLUNI1) ITEXPI = ITSAV c c call myusrout() c RETURN *=== End of subroutine Usrout =========================================* END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: UpdateCurrPart C C This is not a standard FLUKA user routine. It manages the logic to recognize new created particles and to manage C the event history. It is called by the various entries in MGDRAW routine. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c------------------------------------------------------------------------- subroutine UpdateCurrPart(MREG,ICPART,IBPART,AMPART,ICODE, & XX,YY,ZZ,JFLAG) c-------------------------------------------------------------------------- c INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(TRACKR)' INCLUDE '(CASLIM)' INCLUDE '(BEAMCM)' INCLUDE '(RDCYCM)' INCLUDE '(FLKSTK)' INCLUDE '(GENSTK)' INCLUDE '(FHEAVY)' INCLUDE '(PAPROP)' include 'mgdraw.inc' integer JFLAG, JAPART, JZPART logical OLDONE c c JFLAG = 1 -> init tracking, JFLAG = 0 -> tracking of the particle during the event c integer idfluka_now if(jtrack.eq.-1) then return endif c idfluka_now = ISPUSR(MKBMX2) c c ______________ START INIT ________________________ c c Event initialinzing of the kinematics tracking c if(jflag.eq.JSTART) then IONID = ILOFLK (NPFLKA) idfluka(1) = 1 nump = 1 idcurr = 1 idpa(idcurr) = 0 tempo(idcurr) = sngl( agestk(npflka)) tof(idcurr) = 0. trlen(idcurr) = 0. numreg(idcurr) = nrgflk(npflka) igen(idcurr) = 1 idead(idcurr) = 0 idflukamax = 1 if( ionid .GE. 100000 ) then JAPART = MOD ( ILOFLK (NPFLKA), 100000 ) / 100 JZPART = MOD ( ILOFLK (NPFLKA), 10000000 ) / 100000 AAA = JAPART ZZZ = JZPART AMPART = JAPART * AMUC12 & + EMVGEV * EXMSAZ (AAA,ZZZ,.TRUE.,IZDUMM) jpa (idcurr) = -2 else if ( abs(ionid) .ge. 10000 ) then call dcdion ( ionid ) JZPART = ichrge (ionid) JAPART = ibarch (ionid) AMPART = am (ionid) jpa (idcurr) = -2 else if ( ionid .lt. -6) then JZPART = icheav(-ionid) JAPART = ibheav (-ionid) AMPART = amnhea (-ionid) jpa (idcurr) = ionid else JZPART = ichrge (ionid) JAPART = ibarch (ionid) AMPART = am (ionid) jpa (idcurr) = ionid endif c vxi (idcurr) = sngl(xflk (npflka)) vyi (idcurr) = sngl(yflk (npflka)) vzi (idcurr) = sngl(zflk (npflka)) vxf (idcurr) = vxi(idcurr) vyf (idcurr) = vyi(idcurr) vzf (idcurr) = vzi(idcurr) c c when decay the momentum is defined <0 ? c ppp = abs(pmoflk (npflka)) px (idcurr) = sngl(txflk (npflka) * ppp) py (idcurr) = sngl(tyflk (npflka) * ppp) pz (idcurr) = sngl(tzflk (npflka) * ppp) pxf(idcurr) = px(idcurr) pyf(idcurr) = py(idcurr) pzf(idcurr) = pz(idcurr) icha(idcurr) = JZPART iba(idcurr) = JAPART amass(idcurr) = sngl(ampart) c if(idbflg.gt.1 .AND. idbflg.lt.4) then write(*,*)' ' write(*,*)'===== UpdateCurrPart : Ev =',numev,' == INIT = ' write(*,*)' ' write(*,*)' ionid= ',ionid,' npflka= ',npflka, & ' Z= ',jzpart,' A= ',jApart write(*,*)'age= ',tempo(idcurr),' reg= ',numreg(idcurr) write(*,*)'vert= ',xflk (npflka),yflk (npflka),zflk (npflka) write(*,*)'pmod= ',pmoflk (npflka) ppp= abs(pmoflk (npflka)) write(*,*)px(idcurr),py(idcurr),pz(idcurr) write(*,*)'mass= ',ampart,' iloflk= ',iloflk(npflka), & ' jtrack= ',jtrack write(*,*)'idcurr= ',idcurr,'idfluka(idcurr)= ', & idfluka(idcurr),' idfluka_now= ',idfluka_now write(*,*)' ' write(*,*)'===== UpdateCurrPart : END INIT = ' write(*,*)' ' endif return else c _____________________ END INIT ________________________ c c normal Kinematic tracking c if(idbflg.gt.1 .AND. idbflg.lt.4) then write(*,*)' ' write(*,*) & '--------- UpdateCurrPart : Ev =',numev,' ------- ' write(*,*)' ' write(*,*)'JFLAG= ',JFLAG,' Icode = ',icode," nump= ",nump write(*,*)'idcurr= ',idcurr,'idfluka(icurr)= ', & idfluka(idcurr),'idfluka_now= ',idfluka_now write(*,*)'jtrack = ',jtrack,' jpa(idcurr)= ',jpa(idcurr) write(*,*)'PART_CHG = ',icpart,' PART_MASS = ',ampart write(*,*)'ptrack = ',ptrack,'px,y,z = ', & ptrack*cxtrck,ptrack*cytrck,ptrack*cztrck write(*,*)'x(0) = ',xtrack(0),ytrack(0),ztrack(0) write(*,*)'x(1) = ',xtrack(1),ytrack(1),ztrack(1) write(*,*)'atrack = ',atrack,' time = ',tempo(idcurr), % ' tof = ',tof(idcurr) write(*,*)'totlen = ',trlen(idcurr),' ctrack = ',ctrack, % ' cmtrack = ',cmtrck endif c c =========: FOLLOWING THE SAME PARTICLE or an OLD ONE ?: =============== c OLDONE = .false. if (idfluka(idcurr).eq.idfluka_now ) then if (icode.ne.219 .or. icode.eq.106 ) then oldone = .true. else if (idbflg.gt.1 .AND. idbflg.lt.4) & write(*,*) '--------- UpdateCurrPart : Icode 1 =',icode endif c c new charged particle below threshold, giving the energy deposition to father c else if( (ptrack.eq.0).and. (JFLAG.eq.JEND) ) then if (icode.ne.219) then oldone = .true. else if (idbflg.gt.1 .AND. idbflg.lt.4) & write(*,*) '--------- UpdateCurrPart : Icode 2 =',icode endif if(idbflg.gt.0 .AND. idbflg.lt.4) then write(*,*) & 'UpdateCurrentParticle ptrack = 0 searching for father' endif if (numint .GT. maxint) & write(*,*) 'UPDATE: numint>maxint',ncase,numint,maxint do ii = 1,numint if(idbflg.gt.0 .AND. idbflg.lt.4) then write(*,*)'int= ',ii, 'x int= ',xint(ii),yint(ii), & zint(ii),' intpa= ',intpa(ii) endif if( (xx.eq.xint(ii)) .and. (yy.eq.yint(ii)) .and. & (zz.eq.zint(ii)) ) then if(idgflg.gt.0) write(*,*)'father= ',ii idcurr = intpa(ii) return endif end do else do kk = 1,nump if (idfluka(kk).eq.idfluka_now ) then if (icode.ne.219) then oldone = .true. else if (idbflg.gt.1 .AND. idbflg.lt.4) & write(*,*) & '--------- UpdateCurrPart : Icode 3 =',icode endif idcurr = kk endif end do endif c if(oldone) then if(JFLAG.eq.JRUN) then vxf(idcurr) = sngl(xtrack(1)) vyf(idcurr) = sngl(ytrack(1)) vzf(idcurr) = sngl(ztrack(1)) elseif(JFLAG.eq.JEND) then vxf(idcurr) = sngl(xx) vyf(idcurr) = sngl(yy) vzf(idcurr) = sngl(zz) else write(*,*)'UpdateCurrentParticle: ERROR -> wrong JFLAG= ', & JFLAG endif pxf(idcurr) = sngl(ptrack*cxtrck) pyf(idcurr) = sngl(ptrack*cytrck) pzf(idcurr) = sngl(ptrack*cztrck) tof(idcurr) = sngl(atrack) - tempo(idcurr) trlen(idcurr) = sngl(cmtrck) else c c =========: FOLLOWING A NEW PARTICLE : =============== c if(idbflg.gt.0 .AND. idbflg.lt.4) then write(*,*)' ' write(*,*) & 'UPdateCurrentParticle= --------- New part found ----------' write(*,*)' ' write(*,*) write(*,*)'jtrack= ',jtrack,' cha= ',icpart, & ' mass= ',ampart write(*,*)'nump= ',nump,' idfluka_now',idfluka_now write(*,*)'JRUN= ',JRUN,'numint= ',numint write(*,*) 'intcode =',icode write(*,*)'ispusr(mkbmx2)= ',ISPUSR(MKBMX2), & 'ispusr(1)= ',ISPUSR(1) write(*,*)'idcurr= ',idcurr,'idfluka(idcurr)= ', & idfluka(idcurr),' idpa(idcurr) = ',idpa(idcurr) write(*,*)'vi= ',xx,yy,zz write(*,*)'pmod= ',ptrack,' ekin= ',etrack-ampart endif c if(nump.eq.maxnump) then write(*,*)'Error: UpdateCurrPart:' write(*,*)'Max particle number exceeded : ev= ',ncase write(*,*)'nump, maxnump: ',nump,maxnump,jtrack return endif c nump = nump + 1 c idcurr = nump idfluka(idcurr) = idfluka_now tempo(idcurr) = sngl(ATRACK) tof(idcurr) = 0. numreg(idcurr) = mreg icha(idcurr) = icpart iba(idcurr) = ibpart amass(idcurr) = sngl(ampart) jpa(idcurr) = jtrack if(JFLAG.eq.JRUN) then vxi(idcurr) = sngl(xtrack(0)) vyi(idcurr) = sngl(ytrack(0)) vzi(idcurr) = sngl(ztrack(0)) vxf(idcurr) = sngl(xtrack(1)) vyf(idcurr) = sngl(ytrack(1)) vzf(idcurr) = sngl(ztrack(1)) elseif(JFLAG.eq.JEND) then vxi(idcurr) = sngl(xx) vyi(idcurr) = sngl(yy) vzi(idcurr) = sngl(zz) vxf(idcurr) = sngl(xx) vyf(idcurr) = sngl(yy) vzf(idcurr) = sngl(zz) endif px(idcurr) = sngl(ptrack*cxtrck) py(idcurr) = sngl(ptrack*cytrck) pz(idcurr) = sngl(ptrack*cztrck) pxf(idcurr) = px(idcurr) pyf(idcurr) = py(idcurr) pzf(idcurr) = pz(idcurr) trlen(idcurr) = sngl(cmtrck) c c particle latching: father, igen c do ii = 1,numint if(idbflg.gt.0 .AND. idbflg.lt.4) then write(*,*)'int= ',ii, 'x int= ',xint(ii),yint(ii), & zint(ii),' intpa= ',intpa(ii) endif if( (xx.eq.xint(ii)) .and. (yy.eq.yint(ii)) .and. & (zz.eq.zint(ii)) ) then if(idgflg.gt.0) write(*,*)'father= ',ii idpa(idcurr) = intpa(ii) igen(idcurr) = igen(intpa(ii)) + 1 endif end do endif ! end of normal tracking c endif c return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C USER ROUTINE: magfld.f C C This subroutine is activated by the MGNFIELD input card. It returns the magnetic field intensity and direction C for the current position and region. This subroutine is called only if the region where the particle is transported C in that moment has a magnetic flag activated in the ASSIGNMA card. In FOOT simulation, this routine reads a C file containing a magnetic field map and interpolates it at run time. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE MAGFLD ( X, Y, Z, T, BTX, BTY, BTZ, B, NREG, IDISC ) INCLUDE '(DBLPRU)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1988-2020 by Alberto Fasso` & Alfredo Ferrari * * All Rights Reserved. * * * * * * Created in 1988 by Alberto Fasso` * * * * * * Last change on 10-Mar-20 by Alberto Fasso` * * SLAC/USA * * * * Input variables: * * x,y,z = current position * * t = current time * * nreg = current region * * Output variables: * * btx,bty,btz = cosines of the magn. field vector * * B = magnetic field intensity (Tesla) * * idisc = set to 1 if the particle has to be discarded * * * *----------------------------------------------------------------------* * INCLUDE '(BLNKCM)' INCLUDE '(CMEMFL)' INCLUDE '(TRACKR)' INCLUDE '(CSMCRY)' * Statement functions: don't include commons belows this line INCLUDE '(SFEMFL)' include "parameters.inc" * INTEGER NX, NY, NZ, NROW INTEGER XDIM, YDIM, ZDIM DOUBLE PRECISION G2T PARAMETER (G2T=1.D-04) PARAMETER (XDIM=50) PARAMETER (YDIM=50) PARAMETER (ZDIM=500) DOUBLE PRECISION XLAT, YLAT, ZLAT DOUBLE PRECISION XDIFF, YDIFF, ZDIFF DOUBLE PRECISION BXLAT, BYLAT, BZLAT DIMENSION XLAT(XDIM), YLAT(YDIM), ZLAT(ZDIM) DIMENSION BXLAT(XDIM,YDIM,ZDIM), BYLAT(XDIM,YDIM,ZDIM), & BZLAT(XDIM,YDIM,ZDIM) INTEGER I0, J0, K0 INTEGER I1, J1, K1 INTEGER ICOUNT DOUBLE PRECISION B00X,B01X,B10X,B11X,B0X,B1X,BXFIN DOUBLE PRECISION B00Y,B01Y,B10Y,B11Y,B0Y,B1Y,BYFIN DOUBLE PRECISION B00Z,B01Z,B10Z,B11Z,B0Z,B1Z,BZFIN LOGICAL LFIRST DATA LFIRST /.TRUE./ SAVE LFIRST IDISC = 0 IF (LFIRST) THEN CALL OAUXFI(mapname,22,'OLD',IERR) WRITE(*,*) 'Opening Magnetic Map: ',mapname READ (22,*) NROW, NX, NY, NZ IF (XDIM.LT.NX .OR. YDIM.LT.NY .OR. ZDIM.LT.NZ) THEN WRITE(*,*)'magfld.f error: problem in arrays dimension' CALL FLABRT('MAGFLD','Error in arrays dimension') ENDIF DO I=1,NX DO J=1,NY DO K=1,NZ ICOUNT = ICOUNT+1 * Filling XB(I), YB(J), ZB(K) arrays. They are the coordinates * of a point in the space mesh, while BXLAT(I,J,K), * BYLAT(I,J,K), BZLAT(I,J,K) arrays are the B components in that point READ (22,*) XLAT(I), YLAT(J), ZLAT(K), & BXLAT(I,J,K), BYLAT(I,J,K), BZLAT(I,J,K) XLAT(I) = XLAT(I)+MagCenterX YLAT(J) = YLAT(J)+MagCenterY ZLAT(K) = ZLAT(K)+MagCenterZ BXLAT(I,J,K) = BXLAT(I,J,K)*G2T BYLAT(I,J,K) = BYLAT(I,J,K)*G2T BZLAT(I,J,K) = BZLAT(I,J,K)*G2T ENDDO END DO END DO IF (ICOUNT.NE.NROW) THEN WRITE(*,*)'Mag field lattice: reading problem' CALL FLABRT('MAGFLD','Error in map readout') ENDIF WRITE(*,*) 'Magnetic Map reading complete' LFIRST = .FALSE. ENDIF IF ( X.LE.XLAT(1) .OR. X.GT.XLAT(NX) .OR. & Y.LE.YLAT(1) .OR. Y.GT.YLAT(NY) .OR. & Z.LE.ZLAT(1) .OR. Z.GT.ZLAT(NZ) ) THEN C GEOMAGNETIC FIELD AT CNAO BXE = -18.28D-06 BYE = 13.61D-06 BZE = -41.80D-06 B = SQRT(BXE**2 + BYE**2 + BZE**2) BTX = BXE/B BTY = BYE/B BTZ = BZE/B ELSE * Finding the XB coordinates of lattice points which are immediately bfore and after * the current X and save indeces (I0 and I1). * The same is repeated for coordinates YB and ZB DO I=1,NX IF (XLAT(I).GE.X) THEN I0=(I-1) I1=(I) GO TO 30 ENDIF END DO WRITE(*,*) ' Warning: I = NX',I,I1,I0,NREG,X 30 DO J=1,NY IF (YLAT(J).GE.Y) THEN J0=(J-1) J1=(J) GO TO 31 ENDIF END DO WRITE(*,*) ' Warning: J = NY',J,J1,J0,NREG,Y 31 DO K=1,NZ IF (ZLAT(K).GE.Z) THEN K0=(K-1) K1=(K) GO TO 32 ENDIF END DO WRITE(*,*) ' Warning: K = NZ',K,K1,K0,NREG,Z * Reference: * https://en.wikipedia.org/wiki/Trilinear_interpolation * 32 XDIFF=(X-XLAT(I0))/(XLAT(I1)-XLAT(I0)) YDIFF=(Y-YLAT(J0))/(YLAT(J1)-YLAT(J0)) ZDIFF=(Z-ZLAT(K0))/(ZLAT(K1)-ZLAT(K0)) * X component B00X=BXLAT(I0,J0,K0)*(1-XDIFF)+BXLAT(I1,J0,K0)*XDIFF B01X=BXLAT(I0,J0,K1)*(1-XDIFF)+BXLAT(I1,J0,K1)*XDIFF B10X=BXLAT(I0,J1,K0)*(1-XDIFF)+BXLAT(I1,J1,K0)*XDIFF B11X=BXLAT(I0,J1,K1)*(1-XDIFF)+BXLAT(I1,J1,K1)*XDIFF B0X=B00X*(1-YDIFF)+B10X*YDIFF B1X=B01X*(1-YDIFF)+B11X*YDIFF * BXFIN=B0X*(1-ZDIFF)+B1X*(ZDIFF) * Y component B00Y=BYLAT(I0,J0,K0)*(1-XDIFF)+BYLAT(I1,J0,K0)*XDIFF B01Y=BYLAT(I0,J0,K1)*(1-XDIFF)+BYLAT(I1,J0,K1)*XDIFF B10Y=BYLAT(I0,J1,K0)*(1-XDIFF)+BYLAT(I1,J1,K0)*XDIFF B11Y=BYLAT(I0,J1,K1)*(1-XDIFF)+BYLAT(I1,J1,K1)*XDIFF B0Y=B00Y*(1-YDIFF)+B10Y*YDIFF B1Y=B01Y*(1-YDIFF)+B11Y*YDIFF BYFIN=B0Y*(1-ZDIFF)+B1Y*(ZDIFF) * Z component B00Z=BZLAT(I0,J0,K0)*(1-XDIFF)+BZLAT(I1,J0,K0)*XDIFF B01Z=BZLAT(I0,J0,K1)*(1-XDIFF)+BZLAT(I1,J0,K1)*XDIFF B10Z=BZLAT(I0,J1,K0)*(1-XDIFF)+BZLAT(I1,J1,K0)*XDIFF B11Z=BZLAT(I0,J1,K1)*(1-XDIFF)+BZLAT(I1,J1,K1)*XDIFF B0Z=B00Z*(1-YDIFF)+B10Z*YDIFF B1Z=B01Z*(1-YDIFF)+B11Z*YDIFF BZFIN=B0Z*(1-ZDIFF)+B1Z*(ZDIFF) B=SQRT(BXFIN**2+BYFIN**2+BZFIN**2) BTX=BXFIN/B BTY=BYFIN/B BTZ=BZFIN/B ENDIF RETURN *=== End of subroutine Magfld =========================================* END