************************************************************************* * * SUBROUTINE HSECST * * * store secondaries * * * * ISECO=0: no storing of secondaries * * ISECO=1: stores in temporary stack JSTAK (no info in kine branch) * * ISECO=2: stores secondaries only if not made in Shower converter * * ISECO=3: stores all secondaries * * ISECO=4: stores secondaries only if made within 100 cm of target * * * * last modified on: 21/03/2005 by R.Holzmann * * * ************************************************************************* IMPLICIT NONE #include "geant321/gcking.inc" #include "geant321/gckine.inc" #include "geant321/gctrak.inc" #include "geant321/gctmed.inc" #include "geant321/gclist.inc" #include "geant321/gcvolu.inc" #include "geant321/gcflag.inc" #include "user.inc" #include "kinetups.inc" INTEGER NUMED_DRIFTGAS, NUMED_LEAD, NUMED_WIRES COMMON /SHOWMED/ NUMED_DRIFTGAS, NUMED_LEAD, NUMED_WIRES REAL THET, TOFO, UBBU(5), UBUV(5), PVER(4), VER(3), DIST INTEGER KMEC, IKOR, ITYPA, KSEC, I, J, KPAR, NTADD, ITROLD, KPH INTEGER IT, IV, IADR, IADV, ILOOP, IDAL, IPAR, NVER, NUBU, ISECP LOGICAL LSTORE REAL THOP, PDAL(3,2), PAB(2), PEP DATA UBUV /0.,0.,0.,0.,-1./ DATA UBBU /0.,0.,0.,0.,0./ CHARACTER*4 CKCASE CHARACTER*4 CNAMEC(30), CHVOLU(15), CNAMEC1(30) EQUIVALENCE (CHVOLU,NAMES) EQUIVALENCE (CKCASE,KCASE) EQUIVALENCE (CNAMEC,NAMEC) EQUIVALENCE (CNAMEC1,NAMEC1) REAL ACOSN LSTORE = .TRUE. *** Store only when NOT in lead converter or steel plates of SHOWER *** (i.e. trash all showered particles) c IF(NUMED.EQ.NUMED_LEAD .AND. c + (CHVOLU(5).EQ.'S1PB' .OR. CHVOLU(5).EQ.'S2PB') ) c + LSTORE = .FALSE. IF(ISECO(1).EQ.2) THEN IF(CHVOLU(5).EQ.'S1PB' .OR. CHVOLU(5).EQ.'S2PB' .OR. + CHVOLU(5).EQ.'S1SB' .OR. CHVOLU(5).EQ.'S1ST' .OR. + CHVOLU(5).EQ.'S2SB' .OR. CHVOLU(5).EQ.'S2ST') LSTORE=.FALSE. ENDIF IF(ISECO(1).EQ.4) THEN DIST = SQRT(GPOS(1,1)**2+GPOS(2,1)**2+GPOS(3,1)**2) IF(DIST .GT. 100.) LSTORE=.FALSE. ENDIF DO KSEC=1,NGKINE ! loop over secondaries created in this step IF(ISECO(1).EQ.0) THEN ! do not store anything IFLGK(KSEC) = -1 ELSE IF(ISECO(1).EQ.1) THEN ! store in temporary stack only (=default) IFLGK(KSEC) = 0 ENDIF IFLGK(KSEC) = -1 IF(ISECO(1).EQ.2 .AND. LSTORE) THEN ! new trk stored, if NOT in Shower IFLGK(KSEC) = 1 ELSE IF(ISECO(1).EQ.3) THEN ! store always all secondaries as new trks IFLGK(KSEC) = 1 ELSE IF(ISECO(1).EQ.4 .AND. LSTORE) THEN ! store if <100 cm of target IFLGK(KSEC) = 1 ENDIF ENDDO CALL GSKING(0) ! store secondaries according to flags IFLGK ! (those stored in JSTAK only have track nb. of parent) IF(ISECO(1).LE.1 .OR. .NOT.LSTORE) RETURN CALL GFKINE(ITRA,VER,PVER,IPAR,NVER,UBUV,NUBU) ! retrieve current ! track info DO I=1,NMEC ! find out current mechanism KMEC = LMEC(I) IF(KMEC.LE.30) THEN IF(CKCASE.EQ.CNAMEC(KMEC)) GOTO 10 ! this is the current mechanism ENDIF ENDDO IF(CKCASE.EQ.'ANNI' .AND. CNAMEC(LMEC(1)).EQ.'PAIR') KMEC = 31 c special case with gamma->gamma gamma e- process c (pair creation with immediate annihilation of stopped positron) IF(NGKINE.EQ.2) KMEC = 32 ! electron is also stopped 10 CONTINUE DO KSEC=1,NGKINE ! loop over secondaries and save KINE info IF(IFLGK(KSEC).GT.0) THEN ! this track is stored NTRK = NTRK + 1 IF(NTRK.EQ.MAXTRK) THEN ! panic, nb. of tracks exceeds limit WRITE(6, + '('' HSECST: maxtrk='',I6,'' reached: Event aborted!'')') + MAXTRK ENDIF IF(NTRK.GE.MAXTRK) THEN IEOTRI = 1 ! set /GCFLAG/ abort event flag c NEVENT = NEVENT + 1 RETURN ENDIF TRACK(NTRK) = IFLGK(KSEC) ! new track number PREVTRACK(NTRK) = ITRA ! parent track number IF(TRACK(NTRK).NE.NTRK) + WRITE(6,*) ITRA, ' -> ', TRACK(NTRK), NTRK, KMEC, NUMED DO J=1,3 c TVERT(J,NTRK) = LUNIT*VERT(J) ! vertex of parent (KCKINE) TVERT(J,NTRK) = LUNIT*GPOS(J,KSEC) ! vertex of children (KCKING) ENDDO ! (in mm) MECHANISM(NTRK) = 1000*NUMED+KMEC ! save medium number + ! creation mechanism of track PARTID(NTRK) = INT(GKIN(5,KSEC)) ! particle id DO J=1,3 PMOM(J,NTRK) = 1000.0*GKIN(J,KSEC) ! particle momenta in MeV ENDDO GEN1(NTRK) = UBUV(1) ! propagate generator info 1 GEN2(NTRK) = UBUV(2) ! propagate generator info 2 GEN3(NTRK) = UBUV(3) ! propagate generator info WGHTP(NTRK) = UBUV(4) ! i.e. parent weight, user value USERVAL(NTRK) = UBUV(5) IADR = 0 CALL GSKINU(TRACK(NTRK),5,UBUV,IADR) ! attach generator info ENDIF ! to new track ENDDO RETURN ! comment this line to activate Cerenkov saving part c c 1) need to define ISECP in user.inc, range.F and hgeakey.F c 2) need to find a scheme for track numbering => may need another category c IF(ISECP.GT.0) THEN ! save Cerenkov photons too DO KPH=1,NGPHOT NTRK = NTRK + 1 IF(NTRK.EQ.MAXTRK) THEN ! panic, nb. of tracks exceeds limit WRITE(6, + '('' HSECST: maxtrk='',I6,'' reached: Event aborted!'')') + MAXTRK ENDIF IF(NTRK.GE.MAXTRK) THEN IEOTRI = 1 ! set /GCFLAG/ abort event flag c NEVENT = NEVENT + 1 RETURN ENDIF TRACK(NTRK) = NTRK ! new track number PREVTRACK(NTRK) = ITRA ! parent track number DO J=1,3 TVERT(J,NTRK) = LUNIT*XPHOT(J,KPH) ! vertex in mm ENDDO MECHANISM(NTRK) = 1000*NUMED+KMEC ! save medium number + ! creation mechanism of track PARTID(NTRK) = 50 ! particle id DO J=1,3 PMOM(J,NTRK) = 1000.0*XPHOT(J+3,KPH) ! particle momenta in MeV ENDDO GEN1(NTRK) = UBUV(1) ! propagate generator info GEN2(NTRK) = UBUV(2) ! propagate generator info GEN3(NTRK) = UBUV(3) ! propagate generator info WGHTP(NTRK) = UBUV(4) ! i.e. parent weight, etc. USERVAl(NTRK) = UBUV(5) ENDDO ENDIF RETURN ! comment this line to activate Dalitz part DO I=1,NMEC ! loop over mechanisms now KMEC = LMEC(I) IF(KMEC.GT.30) GOTO 12 IF(CKCASE.EQ.CNAMEC(KMEC)) THEN ! this is the current mechanism IF(NGKINE.EQ.3 .AND. CKCASE.EQ.'DCAZ') THEN ! this is a Dalitz ILOOP = 0 DO KSEC=1,NGKINE ! loop over secondaries again ITYPA = INT(GKIN(5,KSEC)) ! particle id IF(ITYPA.EQ.2 .OR. ITYPA.EQ.3) THEN ! e+ or e- ? IF(ILOOP.EQ.0) IDAL = 1 IF(ILOOP.EQ.1) IDAL = 2 DO IKOR=1,3 PDAL(IKOR,IDAL) = GKIN(IKOR,KSEC) ! momenta ENDDO PAB(IDAL) = SQRT(PDAL(1,IDAL)**2 + + PDAL(2,IDAL)**2 + PDAL(3,IDAL)**2) ILOOP = 1 ENDIF ENDDO PEP = PDAL(1,1)*PDAL(1,2) + PDAL(2,1)*PDAL(2,2) + + PDAL(3,1)*PDAL(3,2) IF(PEP.EQ.0. .OR. PAB(1)*PAB(2).EQ.0.) GOTO 11 THOP = ACOSN(PEP/(PAB(1)*PAB(2)))*57.2958 ! e+e- opening angle PRINT *,'Dalitz Theta(1,2)', THOP, (VECT(J),J=1,3) c CALL HCDIR('//PAWC',' ') c CALL HF1(999,THOP,1.) 11 THOP = -1. ENDIF ENDIF 12 continue ENDDO ! end loop over mechanisms RETURN END c c c REAL FUNCTION ACOSN(ARG) c c return acos(arg) and check for validity of argument c IMPLICIT NONE REAL ARG REAL ARG2 ARG2 = ARG IF(ABS(ARG2).GT.1.0) WRITE(6,*) ' *** ERROR in ACOSN: ARG =', ARG IF(ARG2.LT.-1.0) ARG2 = -1.0 IF(ARG2.GT.1.0) ARG2 = 1.0 ACOSN = ACOS(ARG2) RETURN END