C Updated and converted to m1514.f, December 2013 C C MFILL and WRTSUR modified/extended 23-Oct-2013 by NVM C C Updated August 2, 2013 C C To use ROOT geometry, see "ROOT Geometry Mode" section at the very end of this file C with two examples C 06-Jan-2012 and 29-Apr-2013 C C Updated REG1 (tgeo_blmat -> tgeo_blmat1 as a default) and KILLPTCL (IEXPRF) C 20-Nov-2012 SUBROUTINE MARS1514 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! This portion must not be removed from any MARS software and documentation! ! !-------------------------------------------------------------------------------! ! MARS LICENSING AGREEMENT ! ! ! ! The MARS code is being distributed by Nikolai Mokhov of the Accelerator ! ! Physics Center at Fermilab. The code is distributed free of charge to ! ! non-commercial users. MARS was developed in part with funds provided to ! ! Universities Research Association (URA), Inc.under DOE Contract DE-AC02-76 ! ! CH03000 (currently FRA). A re-distribution of the code is not permitted. ! ! The right is reserved to limit distribution, however the United States ! ! Government has a world wide royalty-free license to use the code. ! ! ! ! Each individual who wants to run the code has to register via mokhov@fnal.gov.! ! It is recommended to apply personally for permission in order to gain access ! ! to updated versions of the code for various platforms. To apply for a code, ! ! please copy or download the User's agreement from the official MARS Web site ! ! http://www-ap.fnal.gov/MARS, insert it into an editor and print it with your ! ! institution's letterhead. Sign and date it, and fax it back to Nikolai Mokhov ! ! at 630-840-6039. ! ! ! ! All materials mentioning use of the MARS code system must display the ! ! following reference: ! ! N.V. Mokhov, "The Mars Code System User's Guide", Fermilab-FN-628 (1995); ! ! N.V. Mokhov, S.I. Striganov, "MARS15 Overview", Fermilab-Conf-07/008-AD (2007); ! in Proc. of Hadronic Shower Simulation Workshop, Fermilab, September 2006, ! ! AIP Conf. Proc. 896, pp. 50-60 (2007); ! http://www-ap.fnal.gov/MARS/. ! ! ! ! Neither the name of FRA/Fermilab nor the names of its contributors may be used! ! to endorse or promote products derived from this software without specific ! ! prior written permission. ! ! ! ! WARRANTY DISCLAIMER ! ! ! ! This software is provided by Universities Research Association, Inc ``As Is''.! ! Universities Research Association Inc., expressly disclaims any express or ! ! implied warranties, including, but not limited to, the implied warranties of ! ! merchantability and fitness for a particular purpose. In no event shall ! ! Universities Research Association, Inc., or the Department of Energy be liable! ! for any direct, indirect, incidental, special, exemplary, or consequential ! ! damages (including, but not limited to, procurement of substitute goods or ! ! services; loss of use, data, or profits; or business interruption) however ! ! caused and on any theory of liability, whether in contract, strict liability, ! ! or tort (including negligence or otherwise) arising in any way out of the use ! ! of this software, even if advised of the possibility of such damage. ! ! ! ! LIABILITY DISCLAIMER ! ! ! ! The Licensee takes this software ``As Is'' and except for liability resulting ! ! from any negligent acts or omissions from Universities Research Association, ! ! Inc., the Licensee agrees to hold harmless and defend Universities Research ! ! Association , Inc and the United States Government, Department of Energy from ! ! any and all direct or consequential damages, costs and expenses, including ! ! attorneys fees arising from injury of any type related to the use of this ! ! Software. ! ! ! ! The views and conclusions contained in the software and documentation are ! ! those of the authors and should not be interpreted as representing official ! ! policies, either expressed or implied, of Universities Research Association / ! ! Fermilab. ! ! ! ! Forward any questions to Nikolai Mokhov at mokhov@fnal.gov ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ****************************************************************** * MARS15 * * m1514.f SOURCE FILE * * * * 1. MAIN & USERS - m1514.f * * 2. BLOCK DATA - m15bldt.f,m15bnab.f * * 3. C-SUBROUTINES - m15cstuff.c,m15gui.c,m15mc.c * * 4. INPUT AND OUTPUT - m15in1.f,m15in2.f,m15out.f * * 5. MARSON - m15mareg.f * * 6. TRANSPORT (h,ems,HI) - m15tr.f,m15trems.f * * 7. TRANSPORT (n) - m15trneu.f,m15trneu-mcnp.f * * 8. FIELD AND SYNCH - m15field.f * * 9. DE/DX - m15dedx.f * * 10. GEOMETRY - m15region.f,m15exg.f * * 11. EVENT GENERATOR - m15treem.f,m15eve.f,m15evepi.f,m15evtgen.f * * 12. ELASTIC - m15elast.f * * 13. X-SECTIONS - m15xsec.f * * 14. MORE PHYSICS - m15ph.f * * 15. CEM2007 - m15cem.f * * 16. NEUTRINO - m15neutrino.f * * 17. DEUTERONS - m15deutron.f * * 18. RADIATION - m15rad.f,m15omega.f * * 19. HBOOK BOOKING/FILLING - m15hist1.f,m15hist2.f,m15tuple.f* * 20. UTILITIES, FFREAD - m15util1.f,m15util2.f * * 21. DUMP,MTUPLE,SRCTERM - m15tuple.f,m15srcterm.f * * 22. LAQGSM - m15treem-laqgsm.f,laqgsm1-6.f* ****************************************************************** C C***************************************************************** C........................................................ C MARS15(14) MAIN PROGRAM C AND C USER-SUPPLIED SUBROUTINES C MIXTUR,BEG1,REG1,REG3,FIELD,LEAK,ALIGN,SAGIT,RFCAVT,EDGEUS, C VFAN,MHSETU,MFILL,WRTSUR,TAGPR,TAGGING,KILLPTCL,BLPROCESS AND C MAD-MARS BEAMLINE BUILDER (MMBLB) USER ROUTINES: MARS2BML,BML2MARS ETC. C C DOUBLE PRECISION VERSION C C REVISION: 03-APR-2014 C----------- C PARTICLE ID: C 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 C p n pi+ pi- K+ K- mu+ mu- g e- e+ pbar pi0 d t He3 He4 num nuam nue nuae C 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 C K0L K0S K0 AK0 LAM ALAM SIG+ SIG0 SIG- nbar KSI0 KSI- OM- sib- si0b sib+ AKSI0 ksib+ omb+ C C JJ is PARTICLE ID: = 1-40 for tabulated "elementary" ones C = 1000*Z1 + A1 - Z1 for A1>4 C----------- C AT E 0: DUMP OUTPUT NTIME TIMES, NUPRI=NSTOP/NTIME C NTIME =-1: WRITE TRACKFIRST x,y,z,dcx,dcy for first track C C EDTO() - TOTAL ENERGY DEPOSITION (DENSITY) AFTER SERVN C C STEPEM< MIN(DELTA Z, DELTA R) MUST BE C C NOB>0: STAR DENSITY, FLUENCE, ENERGY DEPOSITION AND PARTICLE SPECTRA C IN NOB SPECIAL REGIONS C HBOOK: FLUENCE IN (CM), ENERGY DEPOSITION IN (GEV/G/CM3) C C EVERYTHING IS NORMALIZED PER 1 INCIDENT PARTICLE C C VOLUMETRIC: C HISTOGRAM ID = NUMH1+NSG+3*(ICL-1)+100*(IHTYP-1)+1000*(NRE-1) C ID=3, 1003, 2003: Star density (stars/cm3) at p > 0.3 GeV/c C ID=9, 1009, 2009: Residual Dose (Rem/hr) at (30days/1day) at AINT C ID=10,1010, 2010: Prompt flux-to-dose based dose equivalent (mSv/hr) at AINT (p/sec) C ID=210, 1210, 2210: Total energy deposition (GeV/g) C C NSURF>0: SURFACE CROSSING ESTIMATOR (UP TO 100 SURFACES): C SURFACE: C ENERGY SPECTRA: ID = NUMH1+NSG+3*(ICL-1)+100*(IHTYP-1)+1000*(NSUR-1) C X-Y ISOFLUENCE: ID = ID1+1000*(NSUR-1) C ID1=411(n), 412(h+/-), 413(gamma), 414(e+mu) C X-Y ISODOSE: ID = ID1+1000*(NSUR-1) C ID1=415 - ENERGY DEPOSITION (mW/g) at AINT (p/s) C ID1=416 - DOSE EQUIVALENT FTD (mSv/hr) at AINT (p/s) C ID1=417 - DOSE EQUIVALENT DIRECT (mSv/hr) at AINT (p/s) C NEUTRINO ENERGY SPECTRA: ID = ID1+1000*(NSUR-1) C ID1=418(numu), 419(anumu), 420(nue), 421(anue) C C NUMH1=1 (DEFAULT) C C IHTYP=5 SPECTRA, FLUENCE AND DOSE EQUIVALENT C IHTYP=6 TIME SPECTRA IN (TMIN-TMAX) INTERVAL (sec) C HISTOGRAMS ARE UNIFORM IN 80 BINS (in nsec!) C C UP TO 100 CYLINDRICAL (ALONG Z-AXIS) OR X-Y (PERPENDICULAR TO C Z-AXIS) SURFACES C C PARTICLE SPECTRA HISTOGRAMING IN NEBIN=80 BINS: 75 LOG + 5 LIN C NEUTRONS: 28 (1.E-12 - 0.0145) + 52 (0.0145 - E0) GEV C OTHERS: 80 (5.E-4 - E0) GEV C C IF(E0<0.0145 GEV) THEN C NEBIN=28 LOG BINS: C NEUTRONS: 28 (1.E-12 - 0.0145) GEV C OTHERS: 28 (5.E-4 - 0.0145) GEV C C NHSPE=0 - dN/dE, DIVIDED BY DEL=DELTA(E) IN (1/CM2/GEV) C NHSPE=1 - E*dN/dE, FOR ALL PARTICLES IN (1/CM2) C C ALL HISTOGRAMS DIVIDED BY VOLUME, SURFACE AREA, OR/AND DELTA_E ! C C DELVOL = 2*PI*R*(RMA-RMI)*(ZMA-ZMI)/NZBIN/NRBIN = DVOL*R C C C User histograms 601 < ID < 700 C XYZ histograms normalized per AINT (p/s) 701 < ID < 1000 C C C ND >0: LOCAL ESTIMATION IN ND DETECTORS FOR L.E.NEUTRONS C C IPHOTO=0 - SAMPLED PHOTOREACTIONS C IPHOTO=1 - FORCED PHOTOREACTIONS C C RNDM, RM48 etc. replaced with GSLRND ! 07/30/13 C Random number call remains the same: C x = rndm(-1.) C both x and rndm are double precision while "-1." is single precision C C--- EVENT GENERATOR CONTROL --- C *********************** C ICEM 4=0 or/and ILAQ(IM)=0: shielding-like applications C E < 3 GeV: CEM; LAQGSM if A<8 or E < 0.02 or not p,n,pi+,pi-,gamma C 3 < E < 5 GeV: above and inclisive Mix-and-Match C E > 5 GeV: inclusive C C ICEM 4=1 or/and ILAQ(IM)=1 (default): majority of applications, C including heavy-ions, nuclides, DPA at E <= 8 GeV C E < 0.3 GeV: CEM; LAQGSM if A<8 or E < 0.02 or not p,n,pi+,pi-,gamma C 0.3 < E < 0.5 GeV: above and LASQGSM Mix-and-Match C 0.5 < E < 8 GeV: LAQGSM C 8 < E < 10 GeV: LAQGSM and inclusive Mix-and-Match C E > 10 GeV: inclusive C C ICEM 4=2 or/and ILAQ(IM)=2: LAQGSM (heavy-ions, nuclides, DPA at E > 8 GeV) C C--------------------------------- C C******************************************************** C CTRL IVIS IEVT IVOL IHIS IDTR C C*** IVIS=0: REGULAR MONTE CARLO RUN (DEFAULT) C*** IVIS=1: MARS GUI C C*** IEVT=0: REGULAR MONTE CARLO RUN (DEFAULT) C*** IEVT>0: MARS EVENT GENERATOR C C CALL EVTGEN(IEVT,ITDNDY,ITDNDX,NEVBIN) C C IEVT=IM (MATERIAL INDEX) C C ITDNDY: 1 - NO DN/DY CALCULATIONS (DEFAULT) C 2 - DN/DY CALCULATIONS C 3 - DN/DMT CALCULATIONS C C ITDNDX: 1 - NO DN/DE OR DN/DX CALCULATIONS C 2 - DN/DE CALCULATIONS (DEFAULT) C 3 - DN/DX CALCULATIONS C 4 - DN/DX CALCULATIONS (DPMJET LIKE) C C NEVBIN - THE NUMBER OF ENERGY BINS AT ITDNDX=2 (DEFAULT: 50) C C*** IVOL=0: REGULAR MONTE CARLO RUN (DEFAULT) C*** IVOL=1: A SHORT MONTE CARLO RUN TO CALCULATE VOLUMES FOR NON-STANDARD REGIONS C CONTAINED IN THE CYLINDER DEFINED BY THE RZVL CARD: C RZVL R1 R2 Z1 Z2 NVTRIAL C DEFAULT: R1=0, R2=REXT, Z1=ZMIN, Z2=ZMAX, NVTRIAL=3000000 C OUTPUT VOLMC.NON IS READY FOR INCLUSION INTO VFAN C C*** IHIS=0: "STANDARD" HISTOGRAM PALETTE WITH 16 COLORS C*** IHIS=1: "B&W" HISTOGRAM PALETTE WITH 8 COLORS C C*** IDTR=0: REGULAR MONTE CARLO RUN (DEFAULT) C*** IDTR=1: DeTra: Solving decay and transmutation equations C*** IDTR=2: DeTra: Processing DeTra output C******************************************************** C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) LOGICAL IND EXTERNAL BLOCK1,BLOCK2,MASNSG,SIBNAB,MAT_DEDX,BNUDOZ,BD2M EXTERNAL CAMERON,CAMER70,BD1,BD2,BD01,BDEJC,MOLLNIX,PKANASA,BARPOA INCLUDE 'biount.inc' INCLUDE 'fnames.inc' INCLUDE 'hist.inc' COMMON & /BLINT1/IO,IBEAM,IZPRJ,IAPRJ,IEDEP,NTIME,NDET,NOB,NHSPE,NF,NFZ, & NTPL,NDM,ILEN,NSURF,NTOFF,IDPSUR(12) & /LOGIND/IND(20) & /BLCTRL/IVIS,IEVT,IVOL,IHIS,IDTR CCH REAL H,EDGLBM COMMON & /BHBOO3/EDGLBL(5),EDGLBM(3,5),IDGLED(5),NEGLED,NHBK PARAMETER (LUNHIST=40) PARAMETER (NH=8000000) ! also in MHBKINIT COMMON/PAWC/H(NH) CCH PARAMETER (ITDNDY=2) PARAMETER (ITDNDX=2) PARAMETER (NEVBIN=50) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C I/O UNITS: C NRCLOS= 1 READ&CLOSE: GEOM.INP (GREAD) C NRCLOS= 1 READ&CLOSE: *.NDT (FREADN) C NRCLOS= 1 READ&CLOSE: MARACT.INP (MARACT) C NRCLOS= 1 WRITE&CLOSE: DUMP (DUMP) C NRCLOS= 1 WRITE&CLOSE: GEOM.DAT (MAIN,GEOMTS) C NWANS= 2 WRITE: ANSYS.ED (SERVN) C NWNEUN= 3 WRITE: NEUTRINO (DECAY,DECMUN) C UNIT= 4 READ: DIP MAP (SUFI) C UNIT= 6 WRITE: SCREEN (INTERMEDIATE) C UNIT= 7 READ: DPM.EVE (BEG1) C UNIT= 8,9,10,11,27,28,29: CEM and LAQGSM C NWMJL= 12 WRITE: EDMJL.GRA (SERVN) C NWPSI= 13 WRITE: PSINEU.GRA (SERVN,NEUOUT) C NWPSI= 14 WRITE: LENEUTRONS (BEGINN,MARSON) C NREAD= 15 READ: MARS.INP (BEGINN) C NWR= 16 WRITE: MARS.OUT (BEGINN,SERVN,NEUOUT,MUOUT) C NWEGH= 17 WRITE: MUON.EGH (RAYMU) C NWEGH= 18 WRITE: TRACK.PLOT (REGION) C NWEGH= 19 WRITE: VERTEX.PLOT(RAYN,RAYMU,EMTRAC,CLENMA,THRESH) C NWEGH= 20 WRITE: MUON.PLOT (RAYMU) C UNIT= 21 READ&CLOSE: XYZHIS.INP (MHISET, MHSETXYZ) C UNIT= 30 WRITE: fort.30 SIGNALS C UNIT= 36 WRITE: fort.36 VFAN VOLUMES AND FIELD MAP TEST PRINTOUT C LUNHIST=40 WRITE: mars.hbook HBOOK C UNIT= 64 WRITE: TRACKFIRST x,y,z,dcx,dcy for first track at NTIME=-1 C UNIT= 65 WRITE: MTUPLE STANDARD GEOMETRY OUTPUT AT IND(1)=T C UNIT= 66 WRITE: MTUPLE-NON NON-STANDARD GEOMETRY OUTPUT C WRITE: surfdet.test INITIALIZATION PRINTOUT FOR SURFACE DETECTORS C UNIT= 79 WRITE: FIRST 30000 SYNCHROTRON PHOTONS GENERATED C UNIT=81-90 WRITE: fort.81 - fort.90 RESERVED FOR SURFACE DETECTOR FILES c----------------------------- write (*,*) ' ' write (*,*) ' Checking MARS.INP for syntax...' call readanalysis(inew) write (*,*) ' ' if (inew.ne.0) then write (*,*) ' STOP: New MARS.INP_corrected was generated' write (*,*) ' Compare it to original MARS.INP,', & ' move to MARS.INP and run again' write (*,*) ' ' stop else write (*,*) ' Syntax is OK' write (*,*) ' ' endif c----------------------------- CALL MPI_USER_START OPEN (UNIT=NREAD,FILE=INPNAM,STATUS='OLD') OPEN (UNIT=NWR, FILE=OUTNAM,STATUS='UNKNOWN') * CALL PREPDPA1 C=== All initialization work (including RNDM) is done in BEGINN === CALL BEGINN C=== First call to rndm, as x=rndm(-1.), is allowed anywhwere below this line === C Note that rndm is actually GSLRND and C that both x and rndm are double precision while "-1." is single precision CALL MPI_USER_NEWSEEDS() IF (IVIS.NE.0) THEN CALL HLIMIT(NH) CALL STTCL(IVIS,IHIS,INPNAM) ELSE IF(IEVT.GT.0) THEN IF(NHBK.GT.0) THEN CALL MHBKINIT(LUNHIST) END IF CALL EVTGEN(IEVT,ITDNDY,ITDNDX,NEVBIN) IF(NHBK.GT.0) THEN CALL MHBKCLOSE(LUNHIST) END IF ELSE IF(.NOT.IND(16)) CALL INITNT() *** CALL INITEMS( 'EMS.TEST' ) CALL INITEGS5 CALL OPENFILES IF(NOB.GT.0.OR.NSURF.GT.0.OR.NHBK.GT.0) THEN CALL MHBKINIT(LUNHIST) END IF CALL MARSON CALL SERVN IF(NOB.GT.0.OR.NSURF.GT.0.OR.NHBK.GT.0) THEN CALL MHBKCLOSE(LUNHIST) END IF CALL MTUPLE END IF CALL RNG_FREE() CALL MPI_USER_FINISH() RETURN END C------------------------------------------------------------ C------------------------------------------------------------ C------------------------------------------------------------ SUBROUTINE REG3 C...................................................... C RE-DEFINES MATERIAL INDICES FOR STANDARD SECTOR C ONLY FOR STANDARD REGIONS DEFINED IN MARS.INP (NOT IN REG1 !!!) C FOR N <= NFZP = MZ*MR*NF C C REVISION: 05-OCT-1999 C...................................................... IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'tally1.inc' *** Examle ************** * * MATIND(15)=2 * DO L=25,44 * MATIND(L)=5 * END DO * ... ************************* RETURN END C------------------------------------------------------- C---------------------------------------------------- SUBROUTINE SUFI C................................................................ C READS MAGNETIC FIELD MAP IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) RETURN END C---------------------------------------------------- C----------------------------------------------- SUBROUTINE ALIGN(VECT,VOUT,NIB,NIM,IFLAG) C................................................ C USED ONLY IF IND(14)=T C C FICTITIOUS SCATTERING: C DISCRETE AT BOUNDARIES BETWEEN REGIONS WITH SPECIAL NUMBERS NIB AND NIM C NIB - PRIOR CROSSING, NIM - AFTER CROSSING C VECT,VOUT: X,Y,Z,DCX,DCY,DCZ,P C C VECT(1)-VECT(6) CAN BE RE-DEFINED TO VOUT(1)-VOUT(6), C NORMALLY AT IND(4)=T C IF SO, VOUT(1)-VOUT(6) MUST BE FILLED AND IFLAG=1 MUST BE RAISED C----- C CREATED: 1994 BY N.MOKHOV C LAST CHANGE: 11-DEC-2001 BY NVM C----- C................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) DIMENSION VECT(7),VOUT(7) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CC 1 DO L=1,6 CC VOUT(L)=VECT(L) CC END DO RETURN END C---------------------------------------------------- SUBROUTINE SAGIT(CHARGE,STEP,VECT,VOUT,NREG,IFLAG) C....................................................... C USED ONLY IF IND(14)=T C C FICTITIOUS SCATTERING: C ON STEP IN REGION NUMBER NREG C VECT,VOUT: X,Y,Z,DCX,DCY,DCZ,P C C VECT(1)-VECT(6) CAN BE RE-DEFINED TO VOUT(1)-VOUT(6), C NORMALLY AT IND(4)=T C IF SO, VOUT(1)-VOUT(6) MUST BE FILLED AND IFLAG=1 MUST BE RAISED C C----- C CREATED: 1994 BY N.MOKHOV C LAST CHANGE: 11-DEC-2001 BY NVM C----- C....................................................... IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) DIMENSION VECT(7),VOUT(7) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CC 1 DO L=1,6 CC VOUT(L)=VECT(L) CC END DO RETURN END C---------------------------------------------------- SUBROUTINE RFCAVT(JJ,STEP,VECT,VOUT,E,TOFF,NIB,NIM,IRF_FLAG) C................................................ C USED ONLY IF IND(14)=T C C RF CAVITY KICK: C DISCRETE AT BOUNDARIES BETWEEN REGIONS WITH SPECIAL NUMBERS NIB AND NIM C NIB - PRIOR CROSSING, NIM - AFTER CROSSING C C INPUT : JJ, STEP, VECT (X,Y,Z,DCX,DCY,DCZ,P), E, TOFF, NIB, NIM C OUTPUT: VOUT (X,Y,Z,DCX,DCY,DCZ,P), E, IFLAG C C VECT(1)-VECT(7) AND E ARE RE-DEFINED TO VOUT(1)-VOUT(7) AND E C IRF_FLAG=1 MUST BE RAISED IF RF IS PRESENT IN THE SYSTEM, C ALLOWING FOR A POSSIBLE ENERGY INCREASE C IRF_FLAG=2 - E_new <0 C C ICTR - ten parameters from UCTR card for arbitrary user control (default: 10*0) C----- C LAST CHANGE: 10-FEB-2010 BY NVM C----- IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'cmasnsg.inc' COMMON & /BLICTR/ICTR(10) DIMENSION VECT(7),VOUT(7) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CC 1 DO L=1,7 CC VOUT(L)=VECT(L) CC END DO IRF_FLAG=0 RETURN END C---------------------------------------------------- SUBROUTINE EDGEUS(U,V) C.................................................... C EDGE-SCATTERING PROBLEM C.................................................... IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END C--------------------------------------------------- SUBROUTINE VFAN(N,V) C........................................... C INITIALIZE VOLUME(N) ARRAY AND FIND V(N) FOR C "DAMAGED" STANDARD REGIONS, EXTENDED-OVERLAPPED, MCNP AND NON-STANDARD REGIONS C ICTR - ten parameters from UCTR card for arbitrary user control (default: 10*0) C C REVISION: 14-OCT-2003 BY NVM C 17-OCT-2003 BY MAK C 04-DEC-2003 BY NVM C 23-AUG-2004 BY NVM C 13-SEP-2004 BY NVM C 17-NOV-2009 BY NVM C 09-NOV-2010 BY NVM C 06-JAN-2012 BY NVM C 22-JAN-2013 BY NVM C REVISION: 21-FEB-2013 BY NVM C........................................... IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) LOGICAL IND COMMON & /LOGIND/IND(20) & /BLCTRL/IVIS,IEVT,IVOL,IHIS,IDTR & /BLICTR/ICTR(10) INCLUDE 'bexgeo.inc' INCLUDE 'blreg1.inc' INCLUDE 'tally1.inc' DATA NENTER/0/ SAVE NENTER C PARAMETER (PI=3.141592653589793227D+00) C PARAMETER (RCOL=1.70D0) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(NENTER.EQ.0) THEN NENTER=1 N1=NFZPEX ! Number of standard+extended regions (used with non-standard geometry) IF(extShapesOverlapped) THEN N1=NFZP ! Number of standard regions (used with extended overlapped geometry) ELSE IF(IND(5).AND.IND(16)) THEN N1=NFZP ! Number of standard regions (used with MCNP geometry) END IF C=== ROOT geometry === IF(IND(19)) THEN c$$$ call tgeo_blvol(VOLUME(N1+1)) call tgeo_blvol(VOLUME(NFZPEX+1)) END IF C===================== C=== MAD-MARS BEAMLINE BUILDER MMBLB === C=== EXAMPLE OF REAL WORK (uncomment and modify * whereever needed) === * IF (IND(13)) THEN * INITINDEX=N1+1 * CALL SET_CURRENT_BEAMLINE(1) * CALL BLVOL(VOLUME(INITINDEX)) * CALL SET_CURRENT_BEAMLINE(2) * CALL BLVOL(VOLUME(INITINDEX)) * END IF C... TUNNEL WALLS AND BEYOND * DO I=1, TUNNEL_NOF_ZONES * CALL TUNNEL_VOL(I, VOLUME(N1+I+TUNNEL_FIRST_ZONE-1)) * END DO C... BEAM ENCLOSURE DONE WITH MC C================================================================== C PUT HERE THE REAL VOLUMES (CM**3) FOR ALL THE NEEDED C STANDARD "DAMAGED", EXTENDED AND NON-STANDARD REGIONS, C I.E. REDEFINE ANY OF THE PRE-DEFINED VOLUMES(L)=0.D0 (EXTENDED AND NON-STANDARD) C OR PRE-CALCULATED STANDARD VOLUMES. C C For example: C C VOLUME(1) =454.D0 ! standard C VOLUME(7) =1524.D0 ! standard C VOLUME(N1+2) =PI*RCOL*RCOL*180.D0 ! extended C VOLUME(N1+365)=2000.D0 ! non-standard VOLUME(N1+1) =1.D0 ! initial, to be redefined in VOLMC.NON: * INCLUDE 'VOLMC.NON' C============================================================ IF(IVOL.EQ.1) THEN WRITE(*,*)' VOLUME MONTE-CARLO IS RUNNING...' ELSE IF(N1.LT.NCELMX.AND.IVIS.EQ.0.AND.IEVT.EQ.0) THEN DO L=N1+1,NCELMX IF(VOLUME(L).GT.0.D0) GO TO 11 END DO WRITE(*,10) 10 FORMAT(/' Region volumes in extended, MCNP and non-stand' & ,'ard geometry sectors are not defined.'/ & ' Use "OPT" instead of "! OPT" in GEOM.INP, if you ' & ,'are in a pure extended geometry'/' mode without', & ' overlapping regions.'// & ' Otherwise, run volume MC, using', & ' "CTRL 3=1" in MARS.INP, to calculate volumes'/ & ' for the above regions (possibly overlapped)', & ' and put VOLMC.NON content'/' into VFAN routine', & ' directly or with INCLUDE statement', & ' Recompile and un normal Monte Carlo after that', & ' with "CTRL 3=0".') STOP 11 CONTINUE END IF END IF END IF V=VOLUME(N) RETURN END C------------------------------------------------- SUBROUTINE MHSETU C........................................................ C SET UP HISTOGRAM ARRAYS C HISTOGRAM ENTRY FOR USER-DEFINED HISTOGRAMING C C HISTOGRAM ID AVAILABLE: 601 < ID < 700 C C HISTOGRAM TYPE: C IHTYP = 1 - COLLISION C IHTYP = 2 - STEP C IHTYP = 3 - ENERGY DEPOSITION C C----- C CREATED: 15-JUN-2000 BY NVM C LAST CHANGE: 09-JUN-2003 BY NVM C----- C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) ************************************************************* * REMEMBER, HBOOK IS A SINGLE PRECISION ENGINE * DON'T FORGET THE 'REAL' DECLARATIONS SUCH AS: * * REAL AA,ELB,X1,X2,Y1,Y2 * CALL HBOOKB(ID,AA,NEB,ELB,0.) * * CALL HBOOK2(ID,TITLE,NX,X1,X2,NY,Y1,Y2,0.) ************************************************************* C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END C------------------------------------------------- SUBROUTINE MFILL(IHTYP,NREG,IM,JJ,E1,E2,DELE,W,X1,Y1,Z1,X2,Y2,Z2, & DCX,DCY,DCZ,STEP,TOF,NI,IDPRC,NREG2) C........................................................ C HISTOGRAM ENTRY FOR USER-DEFINED HISTOGRAMING C C =============================================================== C For neutrinos, this routine is called only at neutrino production C vertex, not at their transport (see Example-7 below). C C AT SUCH A CALL: C C IHTYP=1 C jj=18,19,20 or 21; IDPRC=53, and KORIG=1,16,22 or 23 (see definitions below) C X1=X2, Y1=Y2, Z1=Z2 are the neutrino-production vertex coordinates (STEP=0) C NREG and IM are the number and material index of the region this vertex belongs to C DCX, DCY, DCZ are neutrino direction cosines C W, E1=E2 and TOF are neutrino weight, energy (GeV) and time-of-flight (sec) C NI is a history number C C NEUTRINO PARENT INFO: C C IORIG is parent's ID C KORIG is decay type C EORIG and WORIG are its kinetic energy (GeV) and weight C DCX2, DCY2 AND DCZ2 are parent's direction cosines C C =============================================================== C C HISTOGRAM ID AVAILABLE: 501 < ID < 700 C C CALL TYPE: C IHTYP = 1 - COLLISION C IHTYP = 2 - STEP ("TRACK-LENGTH") C IHTYP = 3 - ENERGY DEPOSITION (LOCAL OR ON THE STEP) C C If there was a boundary crossing, i.e. NREG2.ne.NREG, then C first point (X1,Y1,Z1) belongs to NREG, second point (X2,Y2,Z2) belongs C to NREG2 and is very close to the boundary. C "IHTYP=2" condition must be used to write a file or C Russian Roulette on surface (see examples below). C C PARTICLE ID: C 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 C p n pi+ pi- K+ K- mu+ mu- g e- e+ ap pi0 d t He3 He4 num nuam nue nuae C 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 C K0L K0S K0 AK0 LAM ALAM SIG+ SIG0 SIG- ANEU KSI0 KSI- OM- ASI+ ASI0 ASI- AKSI0 AKSI- AOM- C C JJ is PARTICLE ID: = 1-40 for tabulated "elementary" ones C = 1000*Z1 + A1 - Z1 for A1>4 C C NREG - REGION NUMBER WHERE COLLISION POINT OR STEP BELONG TO C NREG2 - REGION NUMBER FOR THE ADJACENT REGION (in case of boundary crossing) C STEP - STEP (cm), always belongs to region NREG (with material index IM) C C E1 - ENERGY AT STEP START (GeV) C E2 - ENERGY AT STEP END (GeV) C DELE - ENERGY DEPOSITED (GeV), LOCALLY OR ON THE STEP C W - STATISTICAL WEIGHT C X1,Y1,Z1 - COORDINATES AT THE STEP START (cm) C X2,Y2,Z2 - COORDINATES AT THE STEP END (cm) C DCX,DCY,DCZ - DIRECTION COSINES OF THE STEP, i.e. AT ITS BEGINNING C DCX2,DCY2,DCZ2 - DIRECTION COSINES AT THE STEP END (discrete) OF /BLSCTR/ C TOF - TIME-OF-FLIGHT AT THE STEP END (sec) C NI - HISTORY NUMBER C IDPRC - PROCESS ID C P1 and P2 - MOMENTA (GEV/C) AT STEP START AND AND, RESPECTIVELY, OF /BLSCTR/ C C IHTYP=1: E1=E2, DELE=0, STEP=0, X1=X2, Y1=Y2, Z1=Z2 C IHTYP=3: E1#E2, DELE>0; STEP=0 (X1=X2, Y1=Y2, Z1=Z2) OR STEP>0 (X1#X2, Y1#Y2, Z1#Z2) C * ROUTINE IDPRC PROCESS *************************************************************************** * MARSON, 1 DPA * MARSNU 2 RECOIL NUCLEUS LOCAL DEPOSITION (NON-LAQGSM) * 3 Local deposition of heavy ions at AA-vertex (NON-LAQGSM) * 4 STAR DENSITY * 5 FISSION CEM LOCAL (NON-LAQGSM) * 6 FISSION >5 GEV LOCAL (NON-LAQGSM) * 7 HEAVY FRAGMENTS LOCAL DEPOSITION (NON-LAQGSM) * 8 d, t, He3 and He4 LOCAL DEPOSITION (NON-LAQGSM) * EXPECT 9 SURFACE CROSSING BY HADRONS AT IND(6)=T * 10 STAR DENSITY AT IND(6)=T * 11 FLUENCE AT IND(6)=T * THRESH 12 ENERGY DEPOSITION by SUB-THRESHOLD hadrons * 13 SURFACE CROSSING BY NEAR-THRESHOLD HADRONS * 14 FLUENCE BY NEAR-THRESHOLD HADRONS * 15 ENERGY DEPOSITION ON STEP BY NEAR-THRESHOLD HADRONS * RAYN 16 SURFACE CROSSING BY NEUTRAL HADRONS OR IN VACUUM * 17 SURFACE CROSSING BY CHARGED HADRONS IN MATERIAL * 18 ENERGY DEPOSITION ON STEP BY CHARGED HADRONS * 19 FLUENCE BY HADRONS and DPA ON STEP * 20 ENERGY DEPOSITION BY SUB-THRESHOLD e+e- from hadrons * RAYMU 21 e+e- vertex by muon * 22 FLUENCE BY MUONS * 23 SURFACE CROSSING BY MUONS * 24 ENERGY DEPOSITION ON STEP BY MUONS * 25 BREMS VERTEX BY MUONS * 26 NUCL VERTEX BY MUONS * 27 ENERGY DEPOSITION BY SUB-THRESHOLD BREMS OR NUCL PRODUCTS OF MUONS * 28 ENERGY DEPOSITION BY SUB-THRESHOLD MUONS * 29 ENERGY DEPOSITION BY HEAVY FRAGMENTS FROM MUON CAPTURE * 30 MUON DECAY VERTEX * EMTRAC 31 SURFACE CROSSING BY EMS * 32 EMS FLUENCE IN VACUUM OR PHOTON FLUENCE * 33 e+e- FLUENCE AND ENERGY DEPOSITION ON STEP * 34 ENERGY DEPOSITION BY SUB-THRESHOLD EMS * COLLDLT 35 DELTA-ELECTRON VERTEX * trneu-mcnp 36 LE neutron vertex * 37 LE NEUTRON FLUENCE AND DPA ON STEP * 38 LE NEUTRON SURFACE X-ING * 39 LE NEUTRON SURFACE CROSSING AT IND(6)=T * 40 LE NEUTRON FLUENCE AT IND(6)=T * 41 LE NEUTRON LOCAL ED: no new n or gamma generated * 42 LE NEUTRON LOCAL ED: recoil proton below 50 keV (E1=Ep) * 43 LE NEUTRON: n+p -> d+gamma reaction * 44 LE NEUTRON LOCAL ED: vertex non-fission MCNP (E1=Eterm) * 45 LE NEUTRON LOCAL ED: capture on Li or B in non-LAQGSM (E1=Ehi) * 46 LE NEUTRON LOCAL ED: sub-threshold (E1=En) * 47 LE NEUTRON LOCAL ED: SUB-THRESHOLD AT IND(6)=T * 48 LE NEUTRON LOCAL ED: fission * 49 LE NEUTRON LOCAL ED: vertex non-fission BNAB * RAYNUN 50 FLUENCE BY NEUTRINO * THRESH 51 GAS CAPTURE (HE) * COLLIS 52 GAS CAPTURE (LE) * DECAY 53 NEUTRINO PRODUCTION VERTEX * * Origin tagging for hadrons, muons, heavy ions and electromagnetic showers (EMS). * Moved from TAGPR * * Particle JJ of weight W and energy E1 makes a step STEP starting from * point (X,Y,Z) in region NREG. Both this point and the step fully belong * to the region NREG with material index IM. Particle energy at the end * of the step is E2. The current history number is NI. * STEP=0 for sub-threshold particles and neutrinos at production vertex, DCX,DCY,DCZ - neutrino direction * E1=E2 for neutrals and in vacuum. * * This particle (JJ is not 9, 10 or 11) or the EMS (JJ is 9, 10 or 11) was * originated by particle with ID=IORIG of energy EORIG and weight WORIG in * the process KORIG at the point (XORIG,YORIG,ZORIG) in the region NRORIG with * material index IMORIG. * * STEP > 0 for tracking * STEP = 0 for local energy deposition (Delta_E = E1 - E2) * * ICTR - ten parameters from UCTR card for arbitrary user control (default: 10*0) * * IORIG = IKORIG(1,...) - projectile ID * KORIG = IKORIG(2,...) - origin ID * NRORIG= IKORIG(3,...) - origin region number * IMORIG= IKORIG(4,...) - origin material index * * KORIG = 0 - primary beam * 1 - mu/nu, unstable particle decay (pi) * 2 - muons, prompt at hA-vertex * 3 - muons, Bethe-Heitler pair * 4 - muons, e+e- annihilation * 5 - hadrons, hA-vertex * 6 - hadrons, elastic * 7 - hadrons, from muons * 8 - hadrons, unstable particle decay * 9 - hadrons, EMS * 10 - hadrons, recoil LEN * 11 - hadrons, from neutrinos * 12 - EMS, induced by photons from pi0-decay * 13 - EMS, induced by synchrotron photons * 14 - EMS, induced by g,e+,e-, at hA vertex * 15 - EMS, induced by knock-on electrons from muons or hadrons * 16 - neutri, from muon decay * 17 - EMS, induced by prompt e+e- from muons or hadrons * 18 - EMS, induced by brems photons from muon * 19 - EMS, induced by photons from stopped muons * 20 - EMS, induced by photons from low-energy neutrons * 21 - muons, vector mesons * 22 - neutri, K+-, K0, hyperon, He-6 or Ne-18 decays * 23 - neutri, from muon capture C----- C CREATED: 15-JUN-2000 BY NVM C 15-APR-2004 BY NVM C 13-FEB-2006 BY NVM C 03-DEC-2009 BY NVM C 10-JUN-2011 BY NVM C 14-SEP-2011 BY NVM C 26-MAR-2012 BY NVM C 05-DEC-2012 BY NVM C 24-MAY-2013 BY NVM C 23-OCT-2013 BY NVM C 03-APR-2014 BY NVM Comments added on neutrino parent info C----- C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'cmasnsg.inc' INCLUDE 'ptclstack.inc' COMMON & /BPRTAG/XORIG,YORIG,ZORIG,WORIG,EORIG,IORIG,KORIG,NRORIG,IMORIG & /BLICTR/ICTR(10) & /BLDPA/ZRESID,ARESID,ZCNUCL,ACNUCL,TRESID & /BLSCTR/DCX2,DCY2,DCZ2,P1,P2,ALS17,NSG17,IMS17,JJ17 PARAMETER (CLIGHT=29979245800.D0) PARAMETER (AUNIT =0.931494043D0) C MFILL update with EMTH13, 25-Nov-2013 by NVM c$$$ PARAMETER (EMTH13=1.D-13) ************************************************************* Cc Example-1: Tagging-1 C if(ihtyp.ne.3) return C C de=w*dele C if(de.le.0.0d0) return C kkmsk=0 C* D21 C if(nreg.ge.1608.and.nreg.le.1644) then C m=1 C if(nrorig.eq.1603.or.nrorig.eq.1652) then C kkmsk=1 ! up/dw mask as origin C end if C* fill corresponding arrays C ......................... C end if C Cc Example-2: Tagging-2 C C if(nreg.eq.9) then C write(71,100)NI,NREG,IM,JJ,W,E1,E2,STEP, C & XORIG,YORIG,ZORIG,WORIG,EORIG,IORIG,KORIG,NRORIG,IMORIG C 100 format(' NI,NREG,IM,JJ,W,E1,E2,STEP = ',i7,I5,2i3,4(1pe11.4)/ C & ' XORIG,YORIG,ZORIG,WORIG,EORIG,IORIG,KORIG,NRORIG,IMORIG' C & ,'= ',5(1pe11.4),2i3,2i5) C end if *c Example-3: fill in histograms declared in MHSETU * * REMEMBER, HBOOK IS A SINGLE PRECISION ENGINE * DON'T FORGET CONVERSIONS OF THE FOLLOWING TYPE: * * REAL EEH,WWH,XL,YL,WH * EEH=REAL(E1) * WWH=REAL(W) * CALL HFILL(ID,EEH,0.,WWH) * ... * CALL HF2(ID,XL,YL,WWH) * * *c Example-4: Write particles crossed the boundary between regions 2 and 3 * if(ihtyp.eq.2.and.nreg.eq.2.and.nreg2.eq.3) then * write(71,100)ni,JJ,E2,W,X2,Y2,Z2,DCX,DCY,DCZ,step * 100 format(2i7,9(1pe13.5)) * end if * *c Example-5: Write particles crossed a specified boundary * if(ihtyp.eq.2.and.abs(x-5.d0).lt.10.d0) then * if(z1.lt.25.d0.and.z2.gt.25.d0) then * write(71,100)ni,JJ,E2,W,X2,Y2,Z2,DCX,DCY,DCZ,step * 100 format(2i7,9(1pe13.5)) * end if * end if * *c Example-6: collect the recoil energy spectrum produced by nuclear elastic *c collisions of a particle JJ on ZCNUCL > 1 * if(idprc.eq.1.and.zcnucl.gt.1.0d0) then * if(zresid.eq.zcnucl.and.aresid.eq.acnucl) then *c add W*E1 to a corresponding array/histogram ! note here: E1=TRESID * end if * end if * *c Example-7: Write neutrino production vertex info c$$$ if(idprc.eq.53) then c$$$** ctof=clight*tof c$$$ c$$$C=== Neutrino: c$$$ px =e1*dcx c$$$ py =e1*dcy c$$$ pz =e1*dcz c$$$ write(71,101)ni,jj,x1,y1,z1,px,py,pz,e1,w c$$$ 101 format(' Neutrino: ni,jj,x1,y1,z1,px,py,pz,e1,w=', c$$$ & i8,i3,8(1pe14.6)) c$$$ c$$$C=== Parent: c$$$ if(iorig.le.nptcl) then c$$$ smu= pm(iorig) c$$$ else if(iorig.eq.2004) then c$$$ smu=AUNIT*6.015122D0 c$$$ else if(iorig.eq.10008) then c$$$ smu=AUNIT*18.000937D0 c$$$ else c$$$ stop' Stop in MFILL, wrong IORIG' c$$$ end if c$$$ c$$$ p=sqrt(eorig*(eorig+2.d0*smu)) c$$$ px1 =eorig*dcx2 c$$$ py1 =eorig*dcy2 c$$$ pz1 =eorig*dcz2 c$$$ c$$$ write(71,102)iorig,korig,px1,py1,pz1,eorig,worig c$$$ 102 format(' Parent: iorig,korig,px1,py1,pz1,eorig,worig=', c$$$ & 2i4,5(1pe14.6)) c$$$ end if * C Example-8: Write Source Term at Z_SRC * * parameter(z_src=21925.6166759d0) * * if(ihtyp.eq.2.and.z1.lt.z_src.and.z2.ge.z_src) then * IF(E1.GT.EMTH13.AND.E2.GT.EMTH13) THEN * * WRITE(81,102)NI,JJ,E2,W,X2,Y2,Z2,DCX2,DCY2,DCZ2 * 102 FORMAT(I11,I7,2(1PE13.5),2(1PE15.7),1PE19.11,3(1PE17.9)) * * end if * end if ************************************************************* C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END C------------------------------------------------- BLOCK DATA HISTO_DB C........................................................ C USER-DEFINED HISTOGRAMING: C HISTOGRAM ID, TITLES, NORMALIZATION AND BOUNDARIES C C HISTOGRAM ID AVAILABLE: 601 < ID < 700 C C----- C CREATED: 14-AUG-2002 BY NVM C LAST CHANGE: 09-JUN-2003 BY NVM C----- C........................................................ END C------------------------------------------------- C---------------------------------------------------- SUBROUTINE TAGPR(NREG,IM,JJ,W,E1,E2,X,Y,Z,STEP) C.................................................. C C Advise: don't use this routine C Its functions were moved to MFILL C * Origin tagging for hadrons, muons, heavy ions and electromagnetic showers (EMS). * * Particle JJ of weight W and energy E1 makes a step STEP starting from * point (X,Y,Z) in region NREG. Both this point and the step fully belong * to the region NREG with material index IM. Particle energy at the end * of the step is E2. The current history number is NI. * STEP=0 for sub-threshold particles. E1=E2 for neutrals and in vacuum. * * This particle (JJ is not 9, 10 or 11) or the EMS (JJ is 9, 10 or 11) was * originated by particle with ID=IORIG of energy EORIG and weight WORIG in * the process KORIG at the point (XORIG,YORIG,ZORIG) in the region NRORIG with * material index IMORIG. * * STEP > 0 for tracking * STEP = 0 for local energy deposition (Delta_E = E1 - E2) * * ICTR - ten parameters from UCTR card for arbitrary user control (default: 10*0) * * PARTICLE ID (JJ): * 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 * p n pi+ pi- K+ K- mu+ mu- g e- e+ pbar pi0 d t He3 He4 num nuam nue nuae * 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 * K0L K0S K0 AK0 LAM ALAM SIG+ SIG0 SIG- nbar KSI0 KSI- OM- sib- si0b sib+ AKSI0 ksib+ omb+ * * JJ is PARTICLE ID: = 1-40 for tabulated "elementary" ones * = 1000*Z1 + A1 - Z1 for A1>4 * * IORIG = IKORIG(1,...) - projectile ID * KORIG = IKORIG(2,...) - origin ID * NRORIG= IKORIG(3,...) - origin region number * IMORIG= IKORIG(4,...) - origin material index * * KORIG = 0 - primary beam * 1 - mu/nu, unstable particle decay (pi) * 2 - muons, prompt at hA-vertex * 3 - muons, Bethe-Heitler pair * 4 - muons, e+e- annihilation * 5 - hadrons, hA-vertex * 6 - hadrons, elastic * 7 - hadrons, from muons * 8 - hadrons, unstable particle decay * 9 - hadrons, EMS * 10 - hadrons, recoil LEN * 11 - hadrons, from neutrinos * 12 - EMS, induced by photons from pi0-decay * 13 - EMS, induced by synchrotron photons * 14 - EMS, induced by g,e+,e-, at hA vertex * 15 - EMS, induced by knock-on electrons from muons or hadrons * 16 - neutri, from muon decay * 17 - EMS, induced by prompt e+e- from muons or hadrons * 18 - EMS, induced by brems photons from muon * 19 - EMS, induced by photons from stopped muons * 20 - EMS, induced by photons from low-energy neutrons * 21 - muons, vector mesons * 22 - neutri, K+-, K0, hyperon, He-6 or Ne-18 decays * 23 - neutri, from muon capture C----- C CREATED: 11-AUG-2003 BY NVM C 28-APR-2006 BY NVM C 26-APR-2007 BY NVM C 17-NOV-2009 BY NVM C 10-JUN-2011 BY NVM C LAST CHANGE: 25-JAN-2012 BY NVM C----- C.................................................. IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'hist.inc' COMMON & /BPRTAG/XORIG,YORIG,ZORIG,WORIG,EORIG,IORIG,KORIG,NRORIG,IMORIG & /BLICTR/ICTR(10) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C Advise: don't use this routine C Its functions were moved to MFILL C RETURN END C---------------------------------------------------- SUBROUTINE EVENTIN(NST0,NI,IO,WO,EO,AC,ZC,ILAQGSM,ICEM) C.................................................. C Access to a stack of secondaries in every inelastic nuclear C interaction in history NI for projectile IO with weight WO C and kinetic energy EO on nucleus with mass=AC and charge=ZC C for model parameters ILAQGSM and ICEM C C Parameters of (NST-NST0) generated secondaries in the event C in COMMON block of ptclstack.inc C JK(k,i) and EWX(l,i), k=1,3, l=1,10, i=1,NST (nst <= nhamax=10000) C C NST is number of secondary particles in stack C JK(1,i) ID C JK(2,i) generation number C JK(3,i) material index C C EWX(1,i) charge C EWX(2,i) weight C EWX(3,i) kinetic energy (GeV) C EWX(4,i) X C EWX(5,i) Y C EWX(6,i) Z C EWX(7,i) direction cosine DCX C EWX(8,i) direction cosine DCY C EWX(9,i) direction cosine DCZ C C IO and JK(1,i): C is PARTICLE ID = 1-40 for tabulated "elementary" ones C = 1000*Z1 + A1 - Z1 for A1>4 C 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 C p n pi+ pi- K+ K- mu+ mu- g e- e+ ap pi0 d t He3 He4 num nuam nue nuae C 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 C K0L K0S K0 AK0 LAM ALAM SIG+ SIG0 SIG- ANEU KSI0 KSI- OM- ASI+ ASI0 ASI- AKSI0 AKSI- AOM- C C----- C CREATED: 06-JUL-2011 BY NVM C 06-JUL-2011 BY NVM C LAST CHANGE: 13-JAN-2012 BY NVM C----- C.................................................. IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'ptclstack.inc' C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - NSTEV=NST-NST0 ************************************************************* C Example-1: c write(81,101)NI,IO,EO,AC,ZC,ILAQGSM,ICEM,NST c c if(nstev.gt.0) then c do i=NST0+1,NST c jj =jk(1,i) c charge=ewx(1,i) c w =ewx(2,i) c e =ewx(3,i) c smu =ptclmass(jj,charge) c write(81,102)jj,smu,charge,w,e,ewx(7,i),ewx(8,i),ewx(9,i) c end do c end if c 101 format(' NI,IO,EO,AC,ZC,ILAQGSM,ICEM,NST= ',2i7,3(1pe12.3),2i2,i5) c 102 format(i7,1pe14.6,3(1pe12.3),3e15.7) c c c Example-2: c c$$$ c$$$ if(nstev.gt.0) then c$$$ do i=NST0+1,NST c$$$ e =ewx(3,i) c$$$ IERR = ICHKNAN(E) c$$$ IF(IERR.NE.0) go to 1 c$$$ end do c$$$ return c$$$ c$$$ 1 do i=NST0+1,NST c$$$ e =ewx(3,i) c$$$ jj =jk(1,i) c$$$ charge=ewx(1,i) c$$$ w =ewx(2,i) c$$$ smu =ptclmass(jj,charge) c$$$ c$$$ write(*,*)' NI,IO,EO,AC,ILAQGSM,ICEM,NSTEV=', c$$$ & NI,IO,EO,AC,ILAQGSM,ICEM,NSTEV c$$$ write(*,*)' jj,smu,charge,w,e,ewx(9,i)=', c$$$ & jj,smu,charge,w,e,ewx(9,i) c$$$ end do c$$$ stop c$$$ c$$$ end if ************************************************************* RETURN END C------------------------------------------------- SUBROUTINE TAGGING(IM,NREG,WEE) C........................................................ C ENERGY DEPOSITION TAGGING C----- C LAST CHANGE: 08-NOV-2000 BY NVM C----- C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'tally2.inc' C- - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(IM.GT.0) THEN IF(WEE.GT.1.D-50) THEN IF(NUMTAG.GT.0.AND.MTAGG.GT.0) THEN DO N = 1, NUMTAG IF(NREG.EQ.NTAGG(N)) THEN ETAG(INTAG,IETAG,N)=ETAG(INTAG,IETAG,N)+WEE ENDIF ENDDO ENDIF ENDIF ENDIF RETURN END C------------------------------------------------- SUBROUTINE KILLPTCL(JJ,W,E,X,Y,Z,TOFF,NREG,IM) C........................................................ C KILL PARTICLE AT THE BEGINNING OF ITS FIRST STEP (E.G., AT PRODUCTION) C UNDER CERTAIN USER-DEFINED CONDITIONS for any of the input parameter C C ICTR - ten parameters from UCTR card for arbitrary user control (default: 10*0) C C INPUT: JJ,W,E,X,Y,Z,TOFF,NREG,IM C OUTPUT: IEXPRF=-2 to kill C----- C CREATED: 13-APR-2005 BY NVM C 18-NOV-2009 BY NVM C 19-JUN-2012 BY NVM C LAST CHANGE: 20-NOV-2012 BY NVM C----- C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON & /blexprf/iexprf & /BLICTR/ICTR(10) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - *** Example 1: * IF(JJ.EQ.9.AND.NREG.EQ.7.AND.E.GT.7.6D0) THEN * IEXPRF=-2 * END IF *** Example 2: * IF(JJ.EQ.1.AND.Z.GT.150.D0.AND.E.GT.5.D0) THEN * IEXPRF=-2 * END IF *** Example 3: * IF(JJ.EQ.2.AND.W.LT.1.D-30.AND.IM.EQ.3) THEN * IEXPRF=-2 * END IF *** Example 4: * IF(TOFF.GT.1.D-5) THEN * IEXPRF=-2 * END IF RETURN END C------------------------------------------------- BLOCK DATA BLPROCESS C................................................. C IPRCEM(K) = 0 - exclusive C IPRCEM(K) = 1 - inclusive with a probability PRCEM(K) C IPRCEM(0) - global control C C Global bias control is currently done via BIAS card with parameters in C COMMON/BLBIAS/PPIKDEC,PMUPRMT,PMUBEHE,PMUGVM,PMUANN,PPHNUC,PPBAR C C----- C CREATED: 04-OCT-2000 BY NVM C 01-JUN-2001 BY NVM C LAST CHANGE: 10-FEB-2006 BY NVM C----- C................................................. * in process.inc : INTEGER NPRCEM,NPRCHT,NPRCMT,NPRCNT,NPRCHA * in process.inc : INTEGER IPRCEM,IPRCHT,IPRCMT,IPRCNT,IPRCHA * in process.inc : DOUBLE PRECISION PRCEM,PRCHT,PRCMT,PRCNT,PRCHA IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INCLUDE 'process.inc' C- - - - - - - - - - - - - - - - - - - - - - - - - - - - C=== cems === C ... ... ... C 1 2 3 ! set EMS global bias flag: ! 0 - exclusive ! 1 - LPB ! 2 - HALF *** DATA IPRCEM/0,NPRCEM*0/ DATA IPRCEM/1,NPRCEM*0/ DATA PRCEM/NPRCEM*0.D0/ C=== hadron transport === C ... ... ... C 1 2 3 DATA IPRCHT/0,NPRCHT*0/ DATA PRCHT/NPRCHT*0.D0/ C=== muon transport === C ... ... ... C 1 2 3 DATA IPRCMT/0,NPRCMT*0/ DATA PRCMT/NPRCMT*0.D0/ C=== low-energy neutron transport === C ... ... ... C 1 2 3 DATA IPRCNT/0,NPRCNT*0/ DATA PRCNT/NPRCNT*0.D0/ C=== hadron-nucleus vertex === C ... ... ... C 1 2 3 DATA IPRCHA/0,NPRCHA*0/ DATA PRCHA/NPRCHA*0.D0/ END C------------------------------------------------- BLOCK DATA HISTDUMP C........................................................ C VOLUME-SURFACE DETECTOR BINNING (Moved to NOBL and NSUR cards) C C SELECTED DUMP REGION NUMBERS (0 < NBMA < 300) C IN ADDITION TO THE STANDARD GEOMETRY SECTOR C----- C CREATED: 21-NOV-2001 BY NVM C 14-FEB-2002 BY NVM C LAST CHANGE: 11-MAY-2010 BY NVM C----- C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON/BLDUMP/NBU(300),NBMA C=== SELECTED DUMP REGION NUMBERS (0 < NBMA < 300) === DATA NBMA/0/ DATA NBU/300*0/ C For example: C DATA NBMA/6/ C DATA NBU/115,116,117,921,1050,4380, C & 294*0/ C END C------------------------------------------------- ***** MAD-MARS BEAMLINE BUILDER (MMBLB) SECTION ***** C LAST CHANGE: 19-APR-2011 BY NVM ***************************************************** SUBROUTINE MARS2BML(POS,W,BLPOS,BLW) C........................................................ C MARS-TO-BEAMLINE TRANSFORMATION C CALLED IF IND(13)=T before BML2MARS C C----- C CREATED: 2002 BY IT C 22-APR-2003 BY NVM C 19-APR-2011 BY NVM&IT C LAST CHANGE: 27-MAR-2012 BY NVM&IT C----- C........................................................ IMPLICIT NONE DOUBLE PRECISION POS(3), W(3), BLPOS(3), BLW(3) INTEGER ISWITCH,ISWITCH0,IT common/biswitch/iswitch DOUBLE PRECISION MADVEC(3) DOUBLE PRECISION XMAD,YMAD,ZMAD,Z1,Z2 EQUIVALENCE (XMAD,MADVEC(1)),(YMAD,MADVEC(2)),(ZMAD,MADVEC(3)) DOUBLE PRECISION MADDIR(3) DOUBLE PRECISION WXMAD,WYMAD,WZMAD EQUIVALENCE (WXMAD,MADDIR(1)),(WYMAD,MADDIR(2)),(WZMAD,MADDIR(3)) C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C=== EXAMPLE OF MMBLB WORK (1) === * Beamlines are in the positive z-direction * parameter(z1= 650.d0) ! origin of beamline(1), no optics/BLB to the left of this * & ! just standard/non-standard/extended geo in global coordinate system * parameter(z2=5050.d0) ! origin of beamline(2) * C=== Example of work with ROOT BLB (MARBL) === c$$$ IF(IND(19)) THEN c$$$ call mars2mad(POS(1),POS(2),POS(3),XMAD,YMAD,ZMAD) c$$$ call mars2mad(W(1),W(2),W(3),WXMAD,WYMAD,WZMAD) c$$$ c$$$ call tgeo_glob2flat(MADVEC,MADDIR) c$$$ c$$$ call mad2mars(XMAD,YMAD,ZMAD,BLPOS(1),BLPOS(2),BLPOS(3)) c$$$C BLPOS(3)= BLPOS(3)-22751.27475 c$$$ call mad2mars(WXMAD,WYMAD,WZMAD,BLW(1),BLW(2),BLW(3)) c$$$ RETURN c$$$ ENDIF iswitch0 = iswitch iswitch = 0 C=== EXAMPLE OF MMBLB WORK (2) === * * if(POS(3).GT.z1) then * iswitch=1 * else if(POS(3).GT.z2) then * iswitch=2 * endif C================================== C=== Don't touch below!!! === if(iswitch0.eq.-1) return if(iswitch.gt.0) then call set_current_beamline(iswitch) call mars2mad(POS(1),POS(2),POS(3),XMAD,YMAD,ZMAD) call mars2mad(W(1),W(2),W(3),WXMAD,WYMAD,WZMAD) call glob2flat(MADVEC,MADDIR) call mad2mars(XMAD,YMAD,ZMAD,BLPOS(1),BLPOS(2),BLPOS(3)) call mad2mars(WXMAD,WYMAD,WZMAD,BLW(1),BLW(2),BLW(3)) else DO IT=1,3 BLPOS(IT)=POS(IT) BLW(IT)=W(IT) END DO end if C========================================== RETURN END C------------------------------------------------- SUBROUTINE BML2MARS(BML_VEC,BML_DIR,MARS_VEC,MARS_DIR) C........................................................ C BEAMLINE-TO-MARS TRANSFORMATION C CALLED IF IND(13)=T after MARS2BML (which finds iswitch) C C----- C CREATED: 2002 BY IT C 08-NOV-2002 BY NVM C LAST CHANGE: 19-APR-2011 BY NVM&IT C----- C........................................................ IMPLICIT NONE DOUBLE PRECISION BML_VEC(3),BML_DIR(3),MARS_VEC(3),MARS_DIR(3) DOUBLE PRECISION CBUF(3),DBUF(3) INTEGER ISWITCH,IT common/biswitch/iswitch C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(iswitch.gt.0) THEN call set_current_beamline(iswitch) CALL mars2mad(BML_VEC(1),BML_VEC(2),BML_VEC(3), & CBUF(1),CBUF(2),CBUF(3)) CALL mars2mad(BML_DIR(1),BML_DIR(2),BML_DIR(3), & DBUF(1),DBUF(2),DBUF(3)) CALL flat2glob(CBUF,DBUF) CALL mad2mars(CBUF(1),CBUF(2),CBUF(3), & MARS_VEC(1),MARS_VEC(2),MARS_VEC(3)) CALL mad2mars(DBUF(1),DBUF(2),DBUF(3), & MARS_DIR(1),MARS_DIR(2),MARS_DIR(3)) else DO IT=1,3 MARS_VEC(IT)=BML_VEC(IT) MARS_DIR(IT)=BML_DIR(IT) END DO end if RETURN END C-------------------------------------------------------------- SUBROUTINE BLINIT(M_MAX, IMUN) C........................................................ C MMBLB INITIALIZATION C CALLED IF IND(13)=T C C----- C CREATED: 2001 BY KRIOL C MODIFIED: 17-OCT-2003 BY MAK C 04-DEC-2003 BY NVM C LAST CHANGE: 24-AUG-2004 BY MAK&NVM C----- C........................................................ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) INTEGER IMUN(*) LOGICAL MPI_USER_IS_MPI_ACTIVE_RANKS EXTERNAL MPI_USER_IS_MPI_ACTIVE_RANKS include 'blreg1.inc' include 'tally1.inc' C=== EXAMPLE OF REAL WORK (uncomment and modify * whereever needed) === * INCLUDE 'SIMPLE_TUNNEL.INC' * INITZONE= 1 * CALL BEAMLINES_IN_USE(2) ! EXTENSION TO SEVERAL BEAMLINES * CALL SET_CURRENT_BEAMLINE(1) c= Angle of rotation (azimuth) about vertical axis (MARS X-axis), between c= Z axis and projection of the beamline direction vector onto (Z,Y) plane. c= A positive angle forms a right hand screw with MARS X-axis. * BL1_THETA= 0.D0 c= Elevation angle, i.e. the angle between beamline direction and its c= projection onto (Z,Y) plane. If it equals to zero, the beamline remains c= in horizontal plane, positive phi correspond to increasing X. * BL1_PHI = 0.D0 c= Roll angle about the beamline direction vector, i.e. angle between line c= formed by intersection of beamline (x,y)-plane with (Z,Y) plane and c= beamline x axis. A positive angle forms a right hand screw with beamline c= direction vector * BL1_PSI= 0.D0 * x_start = 0.d0 * y_start = 0.d0 * z_start = 0.d0 * call blbuild( 'optics1', * & x_start, y_start, z_start, * & bl1_theta, bl1_phi, bl1_psi, * & InitZone) * The function MPI_USER_IS_MPI_ACTIVE_RANKS is one of interface functions * to the MPI code. It returns .true. if MPI is in use. * The code in the example below makes sure that only the master process * is writing to the file if MPI is working. C... print beamline info into the "BL_ELEMENTS_1" file * IF (MPI_USER_IS_MPI_ACTIVE_RANKS(LOCRANK, MASTERRANK)) THEN ! WITH MPI * IF (LOCRANK .EQ. MASTERRANK) THEN ! THE PROCESS IS MASTER PROCESS * CALL BLPRINT('BL_ELEMENTS_1') * END IF * ELSE ! NO MPI * CALL BLPRINT('BL_ELEMENTS_1') * END IF * write(*,*)'REG1:INFO: use beam Ekin=', * & BeamKineticE, * & ' ;beam mass', * & BeamMass * BeamKineticE = 100.0d0 * BeamMass = 0.93827231d0 * call blsete( BeamKineticE, BeamMass, 5.6780d+03 ) * call blzones( nblnzmax ) * if ( nblnzmax .gt. M_MAX ) then * write(*,*) 'beamline 1: ', * & 'nblnzmax as read from blmatmax ', * & nblnzmax * write(*,*) 'is greater then M_MAX', M_MAX * stop 'STOP' * else if ( nblnzmax .lt. M_MAX ) then * write(*,*) 'beamline 1: ', * & 'nblnzmax as read from blmatmax ', * & nblnzmax * write(*,*) 'is less then M_MAX', M_MAX * endif * if ( nblnzmax .gt. 0 ) then * call blmat(IMUN) * else * write(*,*) 'Number of zones for NuMI beamline <= 0' * stop * endif * call set_current_beamline(2) * bl2_theta = 0.0d0 * bl2_phi = 0.0d0 * bl2_psi = 0.0d0 * x_start = 0.d0 * y_start = -300.d0 * z_start = 3000.d0 * InitZone = 2000 * call blbuild( 'optics2', * & x_start, y_start, z_start, * & bl2_theta, bl2_phi, bl2_psi, * & InitZone) C... print beamline info into the "BL_ELEMENTS_1" file * IF (MPI_USER_IS_MPI_ACTIVE_RANKS(LOCRANK, MASTERRANK)) THEN ! WITH MPI * IF (LOCRANK .EQ. MASTERRANK) THEN ! THE PROCESS IS MASTER PROCESS * CALL BLPRINT('BL_ELEMENTS_2') * END IF * ELSE ! NO MPI * CALL BLPRINT('BL_ELEMENTS_2') * END IF * call blsete( BeamKineticE, BeamMass, 5.6780d+03 ) * call blzones( nblnzmax ) * if ( nblnzmax .gt. M_MAX ) then * write(*,*) 'beamline 1: ', * & 'nblnzmax as read from blmatmax ', * & nblnzmax * write(*,*) 'is greater then M_MAX', M_MAX * stop 'STOP' * else if ( nblnzmax .lt. M_MAX ) then * write(*,*) 'beamline 1: ', * & 'nblnzmax as read from blmatmax ', * & nblnzmax * write(*,*) 'is less then M_MAX', M_MAX * endif * if ( nblnzmax .gt. 0 ) then * call blmat(IMUN) * else * write(*,*) 'Number of zones for NuMI beamline <= 0' * stop * endif * do i=1, st_nof_zones * call simple_tunnel_mat(i, IMUN(i+st_first_zone-1)) * end do c... volume names * call set_current_beamline(1) * call blname(VOLNM(NFZPEX+1)) * call set_current_beamline(2) C* call blname(VOLNM(NFZPEX+1)) RETURN END C------------------------------------------------------------ SUBROUTINE BLGEOINIT(BEAMLINE_ID) C........................................................ C MMBLB ELEMENT REGISTRATION C CALLED IF IND(13)=T C C----- C CREATED: 2001 BY KRIOL C MODIFIED: 17-OCT-2003 BY MAK C LAST CHANGE: 04-DEC-2003 BY NVM C----- C........................................................ IMPLICIT NONE INTEGER BEAMLINE_ID C=== EXAMPLE OF REAL WORK (uncomment and modify * whereever needed) === * INCLUDE 'DRIFT.INC' * INCLUDE 'RBENDEPB.INC' * EXTERNAL * & DRIFT_NAME_FUNC, * & DRIFT_INIT_FUNC, * & DRIFT_MAT_FUNC, * & DRIFT_FIELD_FUNC, * & DRIFT_VOL_FUNC, * & DRIFT_GEO_FUNC, * & DRIFT_ZONENAME_FUNC * IF (BEAMLINE_ID .EQ. 1) THEN * CALL MARS_EL_REGISTER( 1, * & DRIFT_NOF_ZONES, * & DRIFT_NAME_FUNC, * & DRIFT_INIT_FUNC, * & DRIFT_GEO_FUNC, * & DRIFT_MAT_FUNC, * & NO_FIELD, * & DRIFT_VOL_FUNC, * & DRIFT_ZONENAME_FUNC ) * CALL MARS_EL_REGISTER( 2, * & RBENDEPB_NOF_ZONES, * & RBENDEPB_NAME_FUNC, * & RBENDEPB_INIT_FUNC, * & RBENDEPB_GEO_FUNC, * & RBENDEPB_MAT_FUNC, * & RBENDEPB_FIELD_FUNC, * & RBENDEPB_VOL_FUNC, * & RBENDEPB_ZONENAME_FUNC ) * ELSE IF(BEAMLINE_ID .EQ. 2) THEN * CALL MARS_EL_REGISTER( 1, * & DRIFT_NOF_ZONES, * & DRIFT_NAME_FUNC, * & DRIFT_INIT_FUNC, * & DRIFT_GEO_FUNC, * & DRIFT_MAT_FUNC, * & NO_FIELD, * & DRIFT_VOL_FUNC, * & DRIFT_ZONENAME_FUNC ) * CALL MARS_EL_REGISTER( 2, * & RBENDEPB_NOF_ZONES, * & RBENDEPB_NAME_FUNC, * & RBENDEPB_INIT_FUNC, * & RBENDEPB_GEO_FUNC, * & RBENDEPB_MAT_FUNC, * & RBENDEPB_FIELD_FUNC, * & RBENDEPB_VOL_FUNC, * & RBENDEPB_ZONENAME_FUNC ) * END IF RETURN END C------------------------------------------------------------ SUBROUTINE TUNNEL_GEO(XLOCMAD,YLOCMAD,ZLOCMAD,XGLMAD,YGLMAD, & ZGLMAD,NZ,BLINDX,ELINDX,ELBEGINS,ELLENGTH,ANGLE,PREVANGLE, & NEXTANGLE) IMPLICIT NONE INTEGER NZ, BLINDX, ELINDX DOUBLE PRECISION XLOCMAD,YLOCMAD,ZLOCMAD,XGLMAD,YGLMAD,ZGLMAD, & ELBEGINS,ELLENGTH,ANGLE, PREVANGLE, NEXTANGLE RETURN END C------------------------------------------------------------ C C=== ROOT Geometry Mode === C C 1. MARS.INP: Add 19=T to INDX card C C 2. Create a ROOT geometry description in c++ files (see below) C C 3. In GNUmakefile: uncomment #ROOT_ENABLED := YES and comment ROOT_ENABLED := NO; C add the file names to SRCS line as: C SRCS := marsmain.f m1514.f tgeo_init.cc C C 4. Define M_MAX in subroutine REG1 above: C a) for blmat in REG1: M_MAX = 1 + N1 + N2, where N1 is the number of materials C in MATER.INP, and N2 is the number of volumes (shapes/regions/zones) C marked "Red" for "detectors" (scoring regions, sensitive volumes) C in tgeo_init.cc; N1=3, N2=2 in this example, therefore M_MAX=6 in REG1 C C b) for blmat1 in REG1: M_MAX = 1 + N1, where N1 is the total number of regions (nodes) under through node numbering C C 5. Materials description for this simple example is given C in MATER.INP file like: c c NEW Test 01/05/12 c1 'AL' c2 'VAC' c3 'Aluminium honeycomb' 0.05114 5 c 26.98000 13.00000 0.97670 c 14.00700 7.00000 0.01746 c 15.99900 8.00000 0.00552 c 39.94800 18.00000 0.00030 c 1.00794 1.00000 0.00002 cstop C C C 6. A very simple example of ROOT geometry, using the above MATER.INP, C is in a tgeo_init.cc file that must be provided in the same C working directory: C c$$$#include "TROOT.h" c$$$#include "TGeoManager.h" c$$$#include "TGeoMaterial.h" c$$$#include "TGeoTube.h" c$$$#include "TGeoBBox.h" c$$$#include "TGeoMatrix.h" c$$$ c$$$#include c$$$extern "C" { c$$$ c$$$ void tgeo_init_() c$$$ { c$$$ // Substances c$$$ TGeoMedium* GVAC = gGeoManager->GetMedium(0); c$$$ TGeoMedium* AL = gGeoManager->GetMedium("AL"); c$$$ TGeoMedium* VAC = gGeoManager->GetMedium(2); c$$$ TGeoMedium* AlHon = gGeoManager->GetMedium("Aluminium honeycomb"); c$$$ c$$$ // Declare Volumes c$$$ TGeoVolume* World = gGeoManager->MakeTube("world",GVAC,0.,25.,50.); c$$$ c$$$ c$$$ TGeoVolume* box = gGeoManager->MakeBox("box",VAC ,10.,10.,10./2.); c$$$ box->SetLineColor(kRed); c$$$ TGeoVolume* cyl = gGeoManager->MakeTube("cyl",AlHon,0.,10.,10./2.); c$$$ cyl->SetLineColor(kRed); c$$$ c$$$ c$$$ // Positioning volumes c$$$ c$$$ World->AddNode(box,1,new TGeoTranslation("BoxTr",0.,0.,10.)); c$$$ World->AddNode(cyl,1,new TGeoTranslation(10.,0., 20.+10./2.)); c$$$ c$$$ gGeoManager->SetTopVolume(World); c$$$ gGeoManager->CloseGeometry(); c$$$ // Save geometry to file, to have 3D view c$$$ gGeoManager->Export("MyGeom.root"); c$$$ } c$$$} C 7. Another simple example of tgeo_init.cc for a vacuum cylinder (MATER.INP contains just C one material 'VAC') 100cm long and 60cm in radius, that starts at x=y=z=0 in a positive C z-direction and is uniformly subdivided into 10 regions: C c$$$#include c$$$#include c$$$#include c$$$#include c$$$#include c$$$#include "TROOT.h" c$$$#include "TMath.h" c$$$#include "TGeoManager.h" c$$$#include "TGeoMaterial.h" c$$$#include "TGeoTube.h" c$$$#include "TGeoVolume.h" c$$$#include "TGeoMatrix.h" c$$$static const Double_t Zmax=100; // cm, this is a top volume (container) that occupies -ZMAX < z < ZMAX and RGetVolume(); c$$$ if(vol->GetNdaughters() == 0) { c$$$ vol->SetLineColor(kRed); c$$$ } c$$$ } c$$$ return; c$$$ } c$$$ void tgeo_init_() c$$$ { c$$$ c$$$ TGeoMedium* VAC = gGeoManager->GetMedium(1); c$$$ c$$$ TGeoVolume* World=gGeoManager->MakeTube("WORLD",VAC,0.,Rmax,Zmax); c$$$ TGeoVolume* cyl1 = gGeoManager->MakeTube("cyl1",VAC,0.,60.,100./2.); c$$$ TGeoVolume* cyl1Z = cyl1->Divide("cyl1Z",3,10,-100./2.,10.); c$$$ c$$$ SetVolAsSensitive(cyl1); c$$$ // TGeoIterator NodeIter(cyl1); c$$$ // while(TGeoNode* CurrNode = NodeIter()) { c$$$ // TGeoVolume* vol=CurrNode->GetVolume(); c$$$ // if(vol->GetNdaughters() == 0) { c$$$ // vol->SetLineColor(kRed); c$$$ // } c$$$ // } c$$$ c$$$ World->AddNode(cyl1,1,new TGeoTranslation("cyl1Tr",0.,0.,100./2.)); c$$$ c$$$ gGeoManager->SetTopVolume(World); c$$$ gGeoManager->CloseGeometry(); c$$$ gGeoManager->Export("MyGeom.root"); c$$$ } c$$$} C--------------------------------------------------------------------------