module modbacksc implicit none real(8):: user real(4):: userE(2) integer(2):: userI(4) equivalence (user, userE(1)), (user, userI(1)) real(4):: cosOfLargestAngle integer(2):: procIndx, bsflag equivalence (userE(1),cosOfLargestAngle) equivalence (userI(3),procIndx) equivalence (userI(4),bsflag) real(8),parameter::EkchminForLargeAng=1d-2 ! 10 MeV real(8),parameter::EknminForLargeAng=50d-6 ! 50 keV real(8),parameter::backscdircos=0. ! cos(90deg) real(4),allocatable::Elossbs(:,:) integer,parameter::maxindx=6 ! usage of user ! cosOfLargestAngle: to save max scattering angle before the particle ! goes backwoafd ! procIndx: index of the proces caused largest scattering angle ! bsflag: flag to judge a partilce has turned to backward at some point ! (after that it may go forward, but flag (-1) is kept end module modbacksc subroutine epbacksc0 ! to be called from uiaev or cal_uiaev !! use modloft use modbacksc implicit none #include "ZepTrack.h" #include "Zmedia.h" #include "Zzzz.h" allocate ( Elossbs(Ncomp,maxindx) ) ! write(*,*) ! * ' # dE(GeV) Comp Had Brem Knoc anih oth '// ! * ' BS_tot%' ! 20 30.95 11.0 0.0 21.0 17.9 1.4 13.9 65. end subroutine epbacksc1 ! to clear incident track "user" ! to be called from """ uafi1ev **** ! use modbacksc implicit none #include "ZepTrack.h" type(epTrack):: aTrack integer:: icon ! extract one track which should be the only 1 incident call epqsTrack(1, aTrack) ! clear %user userI(:) =0 cosOfLargestAngle = 2.0 aTrack%user = user ! return it to the original position call epretTrack(1, aTrack, icon) Elossbs(:,:) = 0. end subroutine epbacksc1 subroutine epbackscUI(info, loc1, loc2 ) use modbacksc implicit none #include "ZepTrackv.h" #include "Zcode.h" integer,intent(in):: info ! Basically, the particle code ! of the particle which made an interaction. For rare ! particles, they may be given a some unique value ! diff. from acutal code (and branch is made to ! "case default" ) ! Exact code may be obtained by seeing Move.track.p.code. integer,intent(in)::loc1,loc2 ! see above ! To get information of a stacked track, use ! call epgetTrack(i, aTrack, icon) ! where i must be loc1<= i <= loc2. ! record /epTrack/ aTrack is the ouput ! integer icon == 0 ==> i is valid ! icon != 0 ==> i is invalid. ! ! The user can access following variables for the interaction point. ! Move.proc: the interaction type. see each "case" below ! Cn: component number ! Move.track: interacting particle info. e.g Move.track.pos.x etc ! is the position info. in local coordinate of Cn. ! to convert it to world coordinate, ! call epl2w(Cn, Move.track.pos, posw) for position ! call epl2wd(Cn, dir, dirtemp) for direction cos ! Move.track.p.code : code ! Move.track.p.fm.p(1:4): 4 momentum ! Media(MediaNo): media info. E.g Media(MediaNo).name is the media name ! (see Epics/epics/Zmedia.h) ! Media(MediaNo).colA (also colZ, colXs): target of nuclear interaction. ! (A,Z) and cross-section (in mb) ! a(i).p.code etc can be used to know produced particle properties. ! (i=1,n). ! type(epTrack):: aTrack type(epDirec):: dirtemp integer icon, i real(8),external:: kdot_product real(8):: coss ! cos of sacattering angle logical:: ok integer:: indx ! process index. ! write(0,*) ' called with info=',info, loc1,loc2, Move.proc ! write(0,*) 'Cn, code sub charge=',Cn, Move.track.p.code, ! * Move.track.p.subcode,Move.track.p.charge ! write(0,*) 'Energy=',Move.track.p.fm.p(4) select case(info) ! you may freely modify "case" ; ! say, case(kpion:kgnuc) instead of case(kpion:kkaon)... ! etc below case(kphoton) ! gamma interaction ! do something for gamma ineraction ! if you want this routine be called again for this ! event, make info =0 else don't touch it. Then, ! this routine will not be called until next event ! simulation starts. ! Possilbe Move.proc values are ! comp : Compton scattering ! pair : pair creation ! phot : photoelectric effect ! coh : coherent scattering ! photop : photo-production of hadrons ! mpair : magnetic pair creation if( Move%proc == "comp" ) then indx = 1 elseif ( Move%proc == "photop" ) then indx = 2 else indx = 6 ! others endif case(kelec) ! electron interaction ! Move.proc ! brem : Bremstrahlung ! knoc : knockon (Moller or Bhabka) ! anih : positron annihilation ! sync : synchroton emission ! Note-- Cerenkov light emission happens along the ! path of a charged particle and dose not ! call this routine if(Move%proc == "brem" ) then indx = 3 elseif(Move%proc == "knoc" ) then indx = 4 elseif( Move%proc == "anih" ) then indx = 5 else indx = 6 endif case(kmuon) ! muon interaction ! knoc: knock-on ! decay: decay may include negative muon capture ! pair ! brem ! nuci : nuclear interaction indx = 6 case(kpion:kkaon) ! pi, K ! knoc ! decay ! coll collision if( Move%proc == "knoc") then indx = 4 else indx = 2 endif case(knuc) ! p,n ! knoc ! coll : anti-prooton /anti-neutron annihilation included ! To wait until nuclear collision takes place ! if( Move.proc /= 'coll' ) then ! info = 0 ! else ! do i = loc1, loc2 ! print stacked products ! call epgetTrack(i, aTrack, icon) ! write(0,*) i, aTrack.p.code ! enddo ! endif if( Move%proc == "knoc" ) then indx = 4 else indx = 2 endif case(kgnuc) ! heavy int ! coll collision if( Move%proc == "knoc") then indx = 4 else indx = 2 endif ! knoc ! coll : ! case(klight) ! light interaction. This is probably ! nonsense for the moment case default ! others. mainly eta decay ! knoc, coll, decay. indx = 6 end select !!!!!! ! write(0,*) ' index=',indx, ' info=',info !!!!!!! do i = loc1, loc2 call epgetTrack(i, aTrack, icon) user = aTrack%user !!!!! ! write(0,*) ' aTrack.. code=',aTrack%p%code ! write(0,*) ' i=',i, ' bsflag=',bsflag, ! * ' Move.proc=', Move%proc !!!!!!!!! if( bsflag /= -1) then ! check if backward directed here; convert direction to world coord. call epl2wd(aTrack%Cn, aTrack%w, dirtemp) !!!!! ! write(0,*) ' dirtemp', dirtemp%x, dirtemp%y, dirtemp%z !!!!! if( dirtemp%z < backscdircos ) then bsflag = -1 ! backscattered here endif ! cos of scattering angle coss= kdot_product(aTrack%w, Move%Track%w) !!!!!!! ! write(0,*) ' coss=',coss, ' cosOfL.A. ', cosOfLargestAngle !!!!!!! if( aTrack%p%charge /= 0 ) then ok=aTrack%p%fm%p(4)-aTrack%p%mass > EkchminForLargeAng else ok=aTrack%p%fm%p(4)-aTrack%p%mass > EknminForLargeAng endif if( ok .and. coss < cosOfLargestAngle ) then ! memorize largest scattering angle cos ! before begin backscateerred (including ! when backscattered) cosOfLargestAngle = coss procIndx = indx !!!!!! ! write(0,*) ' coss=',coss, ' Procindx=', indx ! write(0,*) ' Move%proc=',Move%proc elseif ( procIndx == 0 ) then ! so far no large angle scattering registered. but now backscaterred ! we must set procIndx ! write(0,*)' coss=',coss, ' indxx=', indx, ' ok=',ok ! write(0,*) ' Move%proc=',Move%proc procIndx = indx !!!!!!!!!! endif aTrack%user = user call epretTrack(i, aTrack, icon) !! else ! ! already backscaterred. don't chage largest angle ! any more. endif enddo end function kdot_product(v1,v2) result(ans) implicit none #include "ZepDirec.h" type(epDirec),intent(in):: v1, v2 real(8):: ans ! scaler product ans = v1%x * v2%x + v1%y * v2%y + v1%z * v2%z end subroutine epbackscdE(cn2, aTrack, dE) use modbacksc implicit none #include "ZepTrack.h" #include "Zmedia.h" #include "Zzzz.h" integer,intent(in):: cn2 ! componet number (as defined in cal_userdE ! (must be called from cal_userdE) ! or simply usual component number depending on the definition. ! ( must be called form userdE) type(epTrack),intent(in):: aTrack real(8),intent(in):: dE ! energy deposit !!!!!!!! ! write(0,*) ' cn2=', cn2, ' dE=',dE !!!!!!!!! user = aTrack%user !!!!!!!! ! write(0,*) ' bsflag=', bsflag, ' procIndx=', procIndx, ! * ' max angl=', cosOfLargestAngle !!!!!! if( bsflag == -1 ) then ! cn = aTrack%cn ! cn2 = DetCN(cn)%num Elossbs(cn2, procIndx) = Elossbs(cn2, procIndx) + dE endif end