c*********************************************************************** c Djangoh user procedure HSUSER * c Creates pre-calculated data file 'djangoh_r.dat' * c Used by DJ_RES function (djangoh_r.f) * c Last update: 14.04.2009 (minor changes) * c * c Additional output (for debug purposes): * c * 'djangoh_r.evts' contain number of events generated * c * 'dj_res_axis.txt' contain axis numbers * c * c Basic parameters (search for lines in code): * c DATA THETAMIN_DEGREE/5.0D0/: Pipe polar angle (in degree) * c NPTSX=200 : number of points on x-axis * c NPTSY=200 : number of points on y-axis * c PARAMETER(NPTS_MAX=1000) : maximal value of NPTSX,NPTSY * c NMODEX=1 : linear/logarythmic x-axis scale * c NMODEY=1 : linear/logarythmic y-axis scale * c XPTMIN=1D-3 : minimal x value * c XPTMAX=1D0 : maximal x value * c YPTMIN=1D-3 : minimal y value * c YPTMAX=1D0 : maximal y value * c FBINWX,FBINWY/0.1D0,0.1D0/ : event counting bin width * c * c Note: z-axis in DJANGOH goes along the initial ELECTRON momentum. * c But the variable definition uses opposite direction of z-axis. * c Therefore all {Pz} variables are involved as {-Pz}. * c * c See webpage for details: http://cern.ch/~makarenko/dj_res/ * c * c Nt: [SIGTOT]=nb * c * c TODO: Labels 'vm_born' --> delete born counter here * C TODO: search & delete all NEVTOT_BORN * c TODO: Remove/modify additional tables * c*********************************************************************** C. Data block (Input parameters: NPTSX,NPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NMODEX,NMODEY) C. RPTSX,RPTSY: data points arrays. C. NPTSX,NPTSY,NPTS_MAX: data table size C. XPTMIN,XPTMAX,YPTMIN,YPTMAX: data table margins C. NMODEX,NMODEY: data grid modes (0: linear, 1: logarythmic) C. NSTATUS: initialization flag (0:no data read, 1:ready, 2:error) C. NCTRLSUM: table control sum C. CS_EVT: Cross section per event C vm_born C. NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN,NEVTS_CORR: number of events tables. Used to estimate errors. BLOCK DATA DJ_RES_DATA1 INTEGER NPTSX,NPTSY,NPTS_MAX,NPTS_MAX2,NMODEX,NMODEY,NCTRLSUM PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,CS_EVT,CS_EVT_BORN INTEGER NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN,NEVTS_CORR COMMON/DJRES_DATA1/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN(NPTS_MAX,NPTS_MAX) & ,NEVTS_CORR(NPTS_MAX,NPTS_MAX),CS_EVT,CS_EVT_BORN C. First init C. Note: NPTSX,NPTSY are intialized at user routine first call DATA NPTSX,NPTSY/0,0/ DATA RPTSX/NPTS_MAX*0D0/ DATA RPTSY/NPTS_MAX*0D0/ DATA NSTATUS,NCTRLSUM/0,0/ DATA XPTMIN,XPTMAX,YPTMIN,YPTMAX/0D0,0D0,0D0,0D0/ DATA NEVTS_BORN/NPTS_MAX2*0/ DATA NEVTS_CORR/NPTS_MAX2*0/ DATA NEVTOT_BORN,NEVTOT_CORR/0,0/ DATA CS_EVT/-1D6/ DATA CS_EVT_BORN/-1D6/ END C. Bins information/margins common block C. The actual bin will be {X*(1-FBINWX),X*(1+FBINWX)} C. ==> Bin width is usually greather than distance between points. BLOCK DATA DJ_RES_BINS1 INTEGER NPTS_MAX PARAMETER(NPTS_MAX=1000) REAL*8 FBINWX,FBINWY,XPT_L,XPT_R,XBIN_W,YPT_L,YPT_R,YBIN_W COMMON/DJRES_BINS/FBINWX,FBINWY,XPT_L(NPTS_MAX),XPT_R(NPTS_MAX), & XBIN_W(NPTS_MAX),YPT_L(NPTS_MAX),YPT_R(NPTS_MAX),YBIN_W(NPTS_MAX) DATA FBINWX,FBINWY/0.03D0,0.03D0/ END C.********************************************************************* SUBROUTINE HSUSER(ICALL,XL,YL,Q2L) IMPLICIT DOUBLE PRECISION (A-H,M,O-Z) C IMPLICIT NONE C. Basic paramter: photon registration angle REAL*8 THETAMIN_DEGREE DATA THETAMIN_DEGREE/5.0D0/ C. Copy of COMMON/DJRES_DATA1/ C vm_born INTEGER NPTSX,NPTSY,NPTS_MAX,NPTS_MAX2,NMODEX,NMODEY,NCTRLSUM PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,CS_EVT,CS_EVT_BORN INTEGER NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN,NEVTS_CORR COMMON/DJRES_DATA1/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN(NPTS_MAX,NPTS_MAX) & ,NEVTS_CORR(NPTS_MAX,NPTS_MAX),CS_EVT,CS_EVT_BORN SAVE/DJRES_DATA1/ C. Current run info. ICURR_TABLE=1: Born generation, ICURR_TABLE=2: Corrections INTEGER ICURR_TABLE C. Copy of COMMON/DJRES_BINS/: bins info common block REAL*8 FBINWX,FBINWY,XPT_L,XPT_R,XBIN_W,YPT_L,YPT_R,YBIN_W COMMON/DJRES_BINS/FBINWX,FBINWY,XPT_L(NPTS_MAX),XPT_R(NPTS_MAX), & XBIN_W(NPTS_MAX),YPT_L(NPTS_MAX),YPT_R(NPTS_MAX),YBIN_W(NPTS_MAX) SAVE/DJRES_BINS/ C. Additional events tables (used for debug reasons only) C. #1: Modified JB, theta=0 C. #2: Modified JB, theta=3 C. #3: Modified JB, theta=10 C. #4: Original JB (same to modified JB at theta=90) C. #5: Hadronic variables C. #6: Leptonic variables REAL*8 SIN_3D,SIN_10D PARAMETER(SIN_3D=0.052335956D0, SIN_10D=0.173648178D0) C. NADDS: number of additional tables, NEVTS_ADD_SIZE: allocation size INTEGER NADDSMAX,NADDS,NEVTS_ADD_SIZE PARAMETER(NADDSMAX=6,NEVTS_ADD_SIZE=NPTS_MAX2*NADDSMAX) INTEGER NEVTS_ADD(NPTS_MAX,NPTS_MAX,NADDSMAX) DATA NEVTS_ADD/NEVTS_ADD_SIZE*0/ DATA NADDS/6/ double precision X_ADD(NADDSMAX),Y_ADD(NADDSMAX),Q2_ADD(NADDSMAX) C. ------------ LOCAL DATA ------------ INTEGER NBXL,NBXR,NBYL,NBYR INTEGER IBINS_EXIST INTEGER NDATFILE, NAXISFILE PARAMETER (NDATFILE=31, NAXISFILE=29) INTEGER I1,I2,I3,ITEMP,ITEMP2,ITEMP3,ITEMP4,IDATOK REAL*8 TEMP1,TEMP2,TEMP3,TEMP4,TEMP5 INTEGER NCALL,NCALL_RC,NCA_LOG,NCURR_PERC DATA NCA_LOG/10/ REAL*8 GSP REAL*8 THETAMINSIN,PI double precision XHAD,YHAD,Q2HAD double precision XHJB,YHJB,Q2JB double precision X_ZEUS,Y_ZEUS,Q2_ZEUS double precision X_JBM,Y_JBM,Q2_JBM C. --------- DEBUG variables --------- C. The nummber or 'out-of-range' events generated INTEGER NFBINFAILS(NADDSMAX+1) C. --------- DJANGOH declarations --------- C. declarations for heracles COMMON /HSUNTS/ LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND COMMON /HSOPTN/ INT2(5),INT3(15),ISAM2(5),ISAM3(15), * IOPLOT,IPRINT,ICUT COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20) COMMON /HSELAB/ SP,EELE,PELE,EPRO,PPRO COMMON /HSCUTS/ XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN COMMON /HSGSW1/ MEI,MEF,MQI,MQF,MEI2,MEF2,MQI2,MQF2,MPRO,MPRO2 COMMON /HSIKP/ S,T,U,SS,TS,US,DKP,DKPS,DKQ,DKQS COMMON /HSLABP/ EH,PH,EQH,PQH,ESH,PSH,COSEH,SINEH PARAMETER (NMXHEP=2000) COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP), & PHEP(5,NMXHEP),VHKK(4,NMXHEP) C ICHNN: event type, born(2) or brems(12) COMMON /HSCHNN/ ICHNN COMMON /HSONLY/ IHSONL C REAL RTIME C. declarations for lepto65 COMMON /DJPASS/ NTOT,NPASS,NFAILL COMMON /DJFAIL/ NFAILI(10) COMMON/LEPTOU/CUT(14),LST(40),PARL(30),XSCH,YSCH,W2SCH,Q2SCH,USCH REAL*4 CUT ,PARL ,XSCH,YSCH,W2SCH,Q2SCH,USCH COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) REAL*4 P ,V COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) REAL*4 PARU ,PARJ INTEGER MSTU,MSTJ C. Number of events from sophia COMMON /SPPASS/ NSOPH,NSPOUT,NFAILP,NSPACC C. Variables from original DJANGOH user procedure C NFAILC: errors counters C IFLCNT: init flavour/errors counters INTEGER IFLCNT,NFAILC DIMENSION IFLCNT(-6:6) DIMENSION NFAILC(0:10) LOGICAL LFIRST DATA LFIRST /.TRUE./ SAVE IFLCNT C. Functions used DOUBLE PRECISION DJ_GET_BORN C. --------- DJANGOH user procedure --------- C. ICALL=1 initialization: before generating events C. =2 analysis:after event have been generated C. =3 final:after all generations IF(LFIRST) THEN C. Note: some variables are already inited LFIRST=.FALSE. C. --------- Init basic parameters --------- C. Parameters: NPTSX,NPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NMODEX,NMODEY C. Init points number NPTSX=200 NPTSY=200 C. Init points scale mode (0: linear, 1: logarythmic) NMODEX=1 NMODEY=1 C. Points scale margins XPTMIN=1D-3 XPTMAX=1D0 YPTMIN=1D-3 YPTMAX=1D0 C. Checks IF(NPTSX.GT.NPTS_MAX)NPTSX=NPTS_MAX IF(NPTSY.GT.NPTS_MAX)NPTSY=NPTS_MAX IF((NMODEX.NE.1).AND.(NMODEX.NE.0))NMODEX=1 IF((NMODEY.NE.1).AND.(NMODEY.NE.0))NMODEY=1 IF((XPTMIN.LE.0D0).OR.(XPTMAX.LE.0D0))NMODEX=0 IF((YPTMIN.LE.0D0).OR.(YPTMAX.LE.0D0))NMODEY=0 C. Data points & bins IF(NMODEX.EQ.1)THEN TEMP1=DLOG(XPTMIN) TEMP2=DLOG(XPTMAX) ENDIF DO 51 I1=1,NPTSX IF(NMODEX.EQ.0)THEN IF (I1.EQ.NPTSX) THEN RPTSX(I1)=XPTMAX ELSE RPTSX(I1)=XPTMIN+(XPTMAX-XPTMIN)*(I1-1)/(NPTSX-1) ENDIF ELSE IF (I1.EQ.1) THEN RPTSX(I1)=XPTMIN ELSEIF (I1.EQ.NPTSX) THEN RPTSX(I1)=XPTMAX ELSE RPTSX(I1)=DEXP(TEMP1+(TEMP2-TEMP1)*(I1-1)/(NPTSX-1)) ENDIF ENDIF XPT_L(I1)=RPTSX(I1)*(1D0-FBINWX) XPT_R(I1)=RPTSX(I1)*(1D0+FBINWX) XBIN_W(I1)=XPT_R(I1)-XPT_L(I1) 51 CONTINUE IF(NMODEY.EQ.1)THEN TEMP1=DLOG(YPTMIN) TEMP2=DLOG(YPTMAX) ENDIF DO 52 I1=1,NPTSY IF(NMODEY.EQ.0)THEN IF (I1.EQ.NPTSY) THEN RPTSY(I1)=YPTMAX ELSE RPTSY(I1)=YPTMIN+(YPTMAX-YPTMIN)*(I1-1)/(NPTSY-1) ENDIF ELSE IF (I1.EQ.1) THEN RPTSY(I1)=YPTMIN ELSEIF (I1.EQ.NPTSY) THEN RPTSY(I1)=YPTMAX ELSE RPTSY(I1)=DEXP(TEMP1+(TEMP2-TEMP1)*(I1-1)/(NPTSY-1)) ENDIF ENDIF YPT_L(I1)=RPTSY(I1)*(1D0-FBINWY) YPT_R(I1)=RPTSY(I1)*(1D0+FBINWY) YBIN_W(I1)=YPT_R(I1)-YPT_L(I1) 52 CONTINUE C. Born or Correction flag. IF ((INT3(7).GT.0).OR.(INT3(8).GT.0).OR.(INT3(9).GT.0)) THEN ICURR_TABLE=2 ELSE ICURR_TABLE=1 ENDIF PI=4D0*DATAN(1D0) THETAMINSIN=DSIN(THETAMIN_DEGREE*PI/180D0) C. Kinematic variables GSP=SP-MPRO2-MEI2 C vm_born C... --------- Load previous events file --------- C. ----- Events file format ----- ver.1, One number per line C. (int) IVER :file format version (currently 1) C. (int) fNDPARS :size of parameters array C. (float) fDPARS(NPARAMS) : array of parameters (cuts and other DJANGOH settings) C. (int) NPTSX,NPTSY,NMODEX,NMODEY C. (float) XPTMIN,XPTMAX,YPTMIN,YPTMAX C. (int) NEVTS_BORN,IEVTS_BORN(NPTSX*NPTSY) C. (int) NEVTS_CORR,IEVTS_CORR(NPTSX*NPTSY) C. (int) fNEVTS_AddN (number of additional events tables) C. (int) fNEVTS_Add(NPTSX*NPTSY*fNEVTS_AddN) (additional data) C. (int) NCTRLSUM: currently NPTSX*NPTSY IDATOK=1 NEVTOT_BORN=0 NEVTOT_CORR=0 OPEN(NDATFILE,FILE='djangoh_r.evts',ERR=78,STATUS='OLD') C. TODO: Implement comment lines in events/results file. READ(NDATFILE,*)IVER IF (IVER.EQ.1) THEN READ(NDATFILE,*) ITEMP C IF(ITEMP.LT.11)THEN IF(ITEMP.LT.10)THEN IDATOK=0 WRITE(LUNOUT,*)'Number of parameters to read: ',ITEMP ENDIF DO 60 I1=1,ITEMP READ(NDATFILE,*) TEMP1 C. Parameters order: FBINWX,FBINWY,ICUT,XMIN,XMAX,YMIN,YMAX,Q2MIN,Q2MAX,WMIN,GSP C. Compare these parameters to current ones IF((I1.EQ.1).AND.(DABS(TEMP1-FBINWX).GT.(FBINWX*1D-5)))THEN WRITE(LUNOUT,*)'FBINWX (old,new): ',TEMP1,FBINWX IDATOK=0 ENDIF IF((I1.EQ.2).AND.(DABS(TEMP1-FBINWY).GT.(FBINWY*1D-5)))THEN WRITE(LUNOUT,*)'FBINWY (old,new): ',TEMP1,FBINWY IDATOK=0 ENDIF IF((I1.EQ.3).AND.(DABS(TEMP1-ICUT).GT.1D-6))THEN WRITE(LUNOUT,*)'ICUT (old,new): ',TEMP1,ICUT IDATOK=0 ENDIF IF((I1.EQ.4).AND.(DABS(TEMP1-XMIN).GT.(XMIN*1D-5)))THEN WRITE(LUNOUT,*)'XMIN (old,new): ',TEMP1,XMIN IDATOK=0 ENDIF IF((I1.EQ.5).AND.(DABS(TEMP1-XMAX).GT.(XMAX*1D-5)))THEN WRITE(LUNOUT,*)'XMAX (old,new): ',TEMP1,XMAX IDATOK=0 ENDIF IF((I1.EQ.6).AND.(DABS(TEMP1-YMIN).GT.(YMIN*1D-5)))THEN WRITE(LUNOUT,*)'YMIN (old,new): ',TEMP1,YMIN IDATOK=0 ENDIF IF((I1.EQ.7).AND.(DABS(TEMP1-YMAX).GT.(YMAX*1D-5)))THEN WRITE(LUNOUT,*)'YMAX (old,new): ',TEMP1,YMAX IDATOK=0 ENDIF IF((I1.EQ.8).AND.(DABS(TEMP1-Q2MIN).GT.(Q2MIN*1D-5)))THEN WRITE(LUNOUT,*)'Q2MIN (old,new): ',TEMP1,Q2MIN IDATOK=0 ENDIF IF((I1.EQ.9).AND.(DABS(TEMP1-Q2MAX).GT.(Q2MAX*1D-5)))THEN WRITE(LUNOUT,*)'Q2MAX (old,new): ',TEMP1,Q2MAX IDATOK=0 ENDIF IF((I1.EQ.10).AND.(DABS(TEMP1-WMIN).GT.(WMIN*1D-5)))THEN WRITE(LUNOUT,*)'WMIN (old,new): ',TEMP1,WMIN IDATOK=0 ENDIF IF((I1.EQ.11).AND.(DABS(TEMP1-GSP).GT.(GSP*1D-5)))THEN WRITE(LUNOUT,*)'GSP (old,new): ',TEMP1,GSP IDATOK=0 ENDIF IF((I1.EQ.12).AND.(TEMP1.GT.0D0))CS_EVT=TEMP1 IF((I1.EQ.13).AND.(TEMP1.GT.0D0))CS_EVT_BORN=TEMP1 60 CONTINUE C. (int) NPTSX,NPTSY,NMODEX,NMODEY READ(NDATFILE,*) ITEMP READ(NDATFILE,*) ITEMP2 READ(NDATFILE,*) ITEMP3 READ(NDATFILE,*) ITEMP4 IF ((ITEMP.NE.NPTSX).OR.(ITEMP2.NE.NPTSY).OR. & (ITEMP3.NE.NMODEX).OR.(ITEMP4.NE.NMODEY)) THEN WRITE(LUNOUT,*)'NPTSX (old,new): ',ITEMP,NPTSX WRITE(LUNOUT,*)'NPTSY (old,new): ',ITEMP2,NPTSY WRITE(LUNOUT,*)'NMODEX (old,new): ',ITEMP3,NMODEX WRITE(LUNOUT,*)'NMODEY (old,new): ',ITEMP4,NMODEY IDATOK=0 ENDIF READ(NDATFILE,*) TEMP1 READ(NDATFILE,*) TEMP2 READ(NDATFILE,*) TEMP3 READ(NDATFILE,*) TEMP4 IF ((TEMP1.NE.XPTMIN).OR.(TEMP2.NE.XPTMAX).OR. & (TEMP3.NE.YPTMIN).OR.(TEMP4.NE.YPTMAX)) THEN WRITE(LUNOUT,*)'XPTMIN (old,new): ',TEMP1,XPTMIN WRITE(LUNOUT,*)'XPTMAX (old,new): ',TEMP2,XPTMAX WRITE(LUNOUT,*)'YPTMIN (old,new): ',TEMP3,YPTMIN WRITE(LUNOUT,*)'YPTMAX (old,new): ',TEMP4,YPTMAX IDATOK=0 ENDIF IF (IDATOK.EQ.1) THEN READ(NDATFILE,*)ITEMP IF(ICURR_TABLE.NE.1)THEN NEVTOT_BORN=ITEMP ELSEIF(ITEMP.GT.0)THEN WRITE(LUNOUT,*)'Born events will be regenerated!' CS_EVT_BORN=-1D6 ENDIF DO 62 I1=1,NPTSX DO 62 I2=1,NPTSY READ(NDATFILE,*)ITEMP IF(NEVTOT_BORN.GT.0)NEVTS_BORN(I1,I2)=ITEMP 62 CONTINUE READ(NDATFILE,*)ITEMP IF(ICURR_TABLE.EQ.1)THEN NEVTOT_CORR=ITEMP ELSEIF(ITEMP.GT.0)THEN WRITE(LUNOUT,*)'RC events will be regenerated!' CS_EVT=-1D6 ENDIF DO 63 I1=1,NPTSX DO 63 I2=1,NPTSY READ(NDATFILE,*)ITEMP IF(NEVTOT_CORR.GT.0)NEVTS_CORR(I1,I2)=ITEMP 63 CONTINUE C. Additional tables load (Note on ICURR_TABLE 1=Born, 2=RC) READ(NDATFILE,*)ITEMP IF(ICURR_TABLE.EQ.1)THEN NADDS=ITEMP ELSEIF(ITEMP.GT.0)THEN WRITE(LUNOUT,*)'Additional tables will be regenerated' ENDIF IF (ITEMP.GT.0) THEN DO 64 I3=1,ITEMP DO 64 I1=1,NPTSX DO 64 I2=1,NPTSY READ(NDATFILE,*)ITEMP2 IF(ICURR_TABLE.EQ.1)NEVTS_ADD(I1,I2,I3)=ITEMP2 64 CONTINUE ENDIF READ(NDATFILE,*) NCTRLSUM IF (NCTRLSUM.NE.(NPTSX*NPTSY)) THEN WRITE(LUNOUT,*)'INVALID control sum read:',NCTRLSUM IDATOK=0 ELSE WRITE(LUNOUT,*)'Control sum read:',NCTRLSUM ENDIF ELSE WRITE(LUNOUT,*)'Previous event file uses DIFFERENT settings' WRITE(LUNOUT,*)'The previous results will be overwritten!' ENDIF ELSE WRITE(LUNOUT,*)'Unexpected version of DAT file' WRITE(LUNOUT,*)'The previous results will be overwritten!' IDATOK=0 ENDIF CLOSE(NDATFILE) C. Truncate read data in case of error IF (IDATOK.EQ.0) THEN WRITE(LUNOUT,*)'Error reading events file.' WRITE(LUNOUT,*)'Previous events are ignored.' CS_EVT=-1D6 CS_EVT_BORN=-1D6 IF(NEVTOT_BORN.GT.0)THEN NEVTOT_BORN=0 DO 71 I1=1,NPTSX DO 71 I2=1,NPTSY 71 NEVTS_BORN(I1,I2)=0 ENDIF IF(NEVTOT_CORR.GT.1)THEN NEVTOT_CORR=0 DO 72 I1=1,NPTSX DO 72 I2=1,NPTSY 72 NEVTS_CORR(I1,I2)=0 ENDIF ENDIF IF(NEVTOT_BORN.GT.0)THEN WRITE(LUNOUT,*)'BORN events data succesfully loaded.' WRITE(LUNOUT,*)'N(events)=',NEVTOT_BORN WRITE(LUNOUT,*)'Cross section per event (Born)=',CS_EVT_BORN WRITE(LUNOUT,*)'Total Born c/s [pb]',CS_EVT_BORN*NEVTOT_BORN ENDIF IF(NEVTOT_CORR.GT.0)THEN WRITE(LUNOUT,*)'RC events data succesfully loaded.' WRITE(LUNOUT,*)'N(events)=',NEVTOT_CORR WRITE(LUNOUT,*)'Cross section per event (RC)=',CS_EVT WRITE(LUNOUT,*)'Total cross section [pb]',CS_EVT*NEVTOT_CORR ENDIF GOTO 79 78 WRITE(LUNOUT,*)'No events file exist yet' 79 WRITE(LUNOUT,*)'Loaded BORN:',NEVTOT_BORN,' RC:',NEVTOT_CORR C. Misc. init NCALL=0 NCALL_RC=0 NCURR_PERC=0 DO 91 I1=-6,6 91 IFLCNT(I1)=0 DO 92 I1=0,10 92 NFAILC(I1)=0 DO 93 I1=1,NADDS+1 93 NFBINFAILS(I1)=0 C. Save axis points (for debugging reason) OPEN(NAXISFILE,FILE='dj_res_axis.txt',STATUS='REPLACE') WRITE(NAXISFILE,*)'============================== X' DO 97 I1=1,NPTSX 97 WRITE(NAXISFILE,*)I1,RPTSX(I1) WRITE(NAXISFILE,*)'============================== Y' DO 98 I1=1,NPTSY 98 WRITE(NAXISFILE,*)I1,RPTSY(I1) CLOSE(NAXISFILE) ENDIF C. --------- User action switch --------- GOTO (100,200,300) ICALL 100 CONTINUE C. The initialization is kept by LFIRST flag RETURN 200 CONTINUE C----------------------------------------------------------------------- C. User's action with generated event NCALL=NCALL+1 C vm_born IF (ICURR_TABLE.EQ.1) THEN NEVTOT_BORN=NEVTOT_BORN+1 ELSE IF (ICURR_TABLE.EQ.2) THEN NEVTOT_CORR=NEVTOT_CORR+1 ENDIF C. Caclulate measurable variables {X_ZEUS,Y_ZEUS} C vm_born IF (ICHNN.EQ.6.OR.ICHNN.EQ.7.OR.ICHNN.EQ.8.OR.ICHNN.EQ.12) THEN C. bremsstrahlung case (ICHNN=12 for CC) NCALL_RC=NCALL_RC+1 C. Note: HEPEVT contains 4 final particles: 1.neutrino, 2.quark, 3.photon and 4.proton remnant C. PHEP(X,1):neutrino C. PHEP(X,2):quark(antiquark) C. PHEP(X,3):photon C. PHEP(X,4):proton remnant C. NHEP==6 (4 finals + initials as ISTHEP()='201') C. Warning: PHEP(5,X) is always ZERO!!! Do not use it! IF(NCALL_RC.LT.NCA_LOG)THEN WRITE(LUNOUT,*)'Bremsstrahlung from /HEPEVT/.' DO 207 I1=1,6 WRITE(LUNOUT,*)'_',I1,'_',IDHEP(I1),ISTHEP(I1),JMOHEP(1,I1) & ,PHEP(4,I1),PHEP(5,I1),PHEP(1,I1),PHEP(2,I1),PHEP(3,I1) 207 CONTINUE WRITE(LUNOUT,*)'Bremsstrahlung from /LUJETS/.' DO 208 I1=1,6 WRITE(LUNOUT,*)I1,P(5,I1),P(4,I1),P(1,I1),P(2,I1),P(3,I1) 208 CONTINUE ENDIF C. Hadronic variables (debug) Q2HAD=-TS XHAD=Q2HAD/(YL*GSP-2D0*(EPRO*PHEP(4,3)+PPRO*PHEP(3,3))) YHAD=Q2HAD/XHAD/GSP C. Jacquet-Blondel variables TEMP1=0D0 TEMP2=0D0 TEMP3=0D0 DO 211 I1=1,NHEP IF((ISTHEP(I1).EQ.1).AND.(I1.NE.1).AND.(I1.NE.3))THEN TEMP1=TEMP1+PHEP(1,I1) TEMP2=TEMP2+PHEP(2,I1) TEMP3=TEMP3+(PHEP(4,I1)+PHEP(3,I1)) ENDIF 211 CONTINUE YJB=TEMP3/2D0/EELE Q2JB=(TEMP1*TEMP1+TEMP2*TEMP2)/(1D0-YJB) XJB=Q2JB/GSP/YJB C. Modified Jacquet-Blondel variables (photon momentum is always added) Y_JBM=(TEMP3+PHEP(4,3)+PHEP(3,3))/2D0/EELE Q2_JBM=((TEMP1+PHEP(1,3))**2+(TEMP2+PHEP(2,3))**2)/(1D0-Y_JBM) X_JBM=Q2_JBM/GSP/Y_JBM C. 'ZEUS' variables (observables: JB + photon momentum at Theta>ThetaPipe) TEMP4=DSQRT(PHEP(1,3)**2+PHEP(2,3)**2)/PHEP(4,3) IF(TEMP4.GT.THETAMINSIN)THEN Y_ZEUS=Y_JBM Q2_ZEUS=Q2_JBM X_ZEUS=X_JBM ELSE Y_ZEUS=YJB Q2_ZEUS=Q2JB X_ZEUS=XJB ENDIF C. Other variables for 'additional tables' C. #1: Modified JB, theta=0 Y_ADD(1)=Y_JBM Q2_ADD(1)=Q2_JBM X_ADD(1)=X_JBM C. #2: Modified JB, theta=3 IF(TEMP4.GT.SIN_3D)THEN Y_ADD(2)=Y_JBM Q2_ADD(2)=Q2_JBM X_ADD(2)=X_JBM ELSE Y_ADD(2)=YJB Q2_ADD(2)=Q2JB X_ADD(2)=XJB ENDIF C. #3: Modified JB, theta=10 IF(TEMP4.GT.SIN_10D)THEN Y_ADD(3)=Y_JBM Q2_ADD(3)=Q2_JBM X_ADD(3)=X_JBM ELSE Y_ADD(3)=YJB Q2_ADD(3)=Q2JB X_ADD(3)=XJB ENDIF C. #4: Original JB (same to modified JB at theta=90) Y_ADD(4)=YJB Q2_ADD(4)=Q2JB X_ADD(4)=XJB C. #5: Hadronic variables Y_ADD(5)=YHAD Q2_ADD(5)=Q2HAD X_ADD(5)=XHAD C. #6: Leptonic variables Y_ADD(6)=YL Q2_ADD(6)=Q2L X_ADD(6)=XL C. Log different: XL,YL(leptonic),XHAD,YHAD(hadronic),XSCH,YSCH(scaled),XJB,YJB,X_ZEUS,Y_ZEUS IF(NCALL_RC.LT.NCA_LOG)THEN WRITE(LUNOUT,*)'Bremsstrahlung variables list (x,y,Q2).' WRITE(LUNOUT,*)'Leptonic:',XL,YL,Q2L WRITE(LUNOUT,*)'Hadronic:',XHAD,YHAD,Q2HAD WRITE(LUNOUT,*)'Scaled:',XSCH,YSCH WRITE(LUNOUT,*)'Jacuet-Blondel:',XJB,YJB,Q2JB WRITE(LUNOUT,*)'ZEUS:',X_ZEUS,Y_ZEUS,Q2_ZEUS C WRITE(LUNOUT,*)'------- additional x,y,Q2 ---------' C WRITE(LUNOUT,*)'1:',X_ADD(1),Y_ADD(1),Q2_ADD(1) C WRITE(LUNOUT,*)'2:',X_ADD(2),Y_ADD(2),Q2_ADD(2) C WRITE(LUNOUT,*)'3:',X_ADD(3),Y_ADD(3),Q2_ADD(3) C WRITE(LUNOUT,*)'4:',X_ADD(4),Y_ADD(4),Q2_ADD(4) C WRITE(LUNOUT,*)'5:',X_ADD(5),Y_ADD(5),Q2_ADD(5) C WRITE(LUNOUT,*)'6:',X_ADD(6),Y_ADD(6),Q2_ADD(6) C WRITE(LUNOUT,*)'-----------------------------------' ENDIF ELSE C vm_born C. born case (ICHNN=2 for CC) Y_ZEUS=YL Q2_ZEUS=Q2L X_ZEUS=XL DO 231 I1=1,NADDS Y_ADD(I1)=YL Q2_ADD(I1)=Q2L X_ADD(I1)=XL 231 CONTINUE IF(NCALL.LT.NCA_LOG)THEN C. DEBUG: Test Jacquet-Blondel variables in elastic case TEMP1=0D0 TEMP2=0D0 TEMP3=0D0 DO 232 I1=1,NHEP IF((ISTHEP(I1).EQ.1).AND.(I1.NE.1))THEN TEMP1=TEMP1+PHEP(1,I1) TEMP2=TEMP2+PHEP(2,I1) TEMP3=TEMP3+(PHEP(4,I1)+PHEP(3,I1)) WRITE(LUNOUT,*)I1,':(E,Pz)=',PHEP(4,I1),PHEP(3,I1) WRITE(LUNOUT,*)I1,':(E-Pz)=',PHEP(4,I1)+PHEP(3,I1) ENDIF 232 CONTINUE YJB=TEMP3/2D0/EELE Q2JB=(TEMP1*TEMP1+TEMP2*TEMP2)/(1D0-YJB) XJB=Q2JB/GSP/YJB WRITE(LUNOUT,*)' =========== elastic vars =============' WRITE(LUNOUT,*)'Jacuet-Blondel:',XJB,YJB,Q2JB WRITE(LUNOUT,*)'Elastic variables:',XL,YL,Q2L WRITE(LUNOUT,*)'EELE=',EELE WRITE(LUNOUT,*)'GSP=',GSP ENDIF ENDIF C. --------- Fill the nearest bins --------- C. Find the most left/right bins for {X_ZEUS,Y_ZEUS} IBINS_EXIST=IDJR_FINDAFFBN(X_ZEUS,Y_ZEUS,NBXL,NBXR,NBYL,NBYR) C. Fill the bins IF(IBINS_EXIST.EQ.1)THEN DO 241 I1=NBXL,NBXR DO 241 I2=NBYL,NBYR IF (ICURR_TABLE.EQ.1) THEN NEVTS_BORN(I1,I2)=NEVTS_BORN(I1,I2)+1 ELSEIF (ICURR_TABLE.EQ.2) THEN NEVTS_CORR(I1,I2)=NEVTS_CORR(I1,I2)+1 ENDIF 241 CONTINUE ELSE NFBINFAILS(NADDS+1)=NFBINFAILS(NADDS+1)+1 ENDIF C. Fill additional tables in RC run (ICURR_TABLE==2) IF ((ICURR_TABLE.EQ.2).AND.(NADDS.GE.1)) THEN DO 252 I3=1,NADDS IBINS_EXIST=IDJR_FINDAFFBN(X_ADD(I3),Y_ADD(I3), & NBXL,NBXR,NBYL,NBYR) IF (IBINS_EXIST.EQ.1) THEN DO 251 I1=NBXL,NBXR DO 251 I2=NBYL,NBYR 251 NEVTS_ADD(I1,I2,I3)=NEVTS_ADD(I1,I2,I3)+1 ELSE NFBINFAILS(I3)=NFBINFAILS(I3)+1 ENDIF 252 CONTINUE ENDIF C...count flavors (part of original user procedure) IF (NSPACC.EQ.0) THEN DO 271 I1=-6,6 271 IF (LST(25).EQ.I1) IFLCNT(I1)=IFLCNT(I1)+1 ENDIF C...count errors in DJGEVT (part of original user procedure) DO 272 I1=0,10 272 IF (NFAILI(I1).NE.0) NFAILC(I1)=NFAILC(I1)+1 C. Simple progress watch IF ((DFLOAT(NCALL)/(NEVENT/100)).GE.NCURR_PERC) THEN NCURR_PERC=NCURR_PERC+1 WRITE(LUNOUT,*) 'Progress: ',NCURR_PERC,'%' C. TODO: Flush buffer here ENDIF RETURN 300 CONTINUE C----------------------------------------------------------------------- C. Cross section per event (in pb), note: [SIGTOT]=nb WRITE(LUNOUT,*)'#Events tried: ',NEVENT WRITE(LUNOUT,*)'#User procedure calls: ',NCALL IF (ICURR_TABLE.EQ.1) THEN CS_EVT_BORN=SIGTOT/DFLOAT(NEVENT)*1D3 WRITE(LUNOUT,*)'#Born c/s (pb): ',SIGTOT*1D3 WRITE(LUNOUT,*)'#Born c/s generated: ',CS_EVT_BORN*NEVTOT_BORN IF((NEVTOT_CORR.GT.0).AND.(CS_EVT.GT.0D0))THEN WRITE(LUNOUT,*)'#Total c/s loaded: ',CS_EVT*NEVTOT_CORR ENDIF ELSEIF (ICURR_TABLE.EQ.2) THEN CS_EVT=SIGTOT/DFLOAT(NEVENT)*1D3 WRITE(LUNOUT,*)'#Total c/s (pb): ',SIGTOT*1D3 WRITE(LUNOUT,*)'#Total c/s generated: ',CS_EVT*NEVTOT_CORR IF((NEVTOT_BORN.GT.0).AND.(CS_EVT_BORN.GT.0D0))THEN WRITE(LUNOUT,*)'#Born c/s loaded: ',CS_EVT_BORN*NEVTOT_BORN ENDIF ENDIF IF((NEVTOT_BORN.GT.0).AND.(CS_EVT_BORN.GT.0D0))THEN WRITE(LUNOUT,*)'#Born c/s per event: ',CS_EVT_BORN WRITE(LUNOUT,*)'#Born events: ',NEVTOT_BORN ENDIF IF((NEVTOT_CORR.GT.0).AND.(CS_EVT.GT.0D0))THEN WRITE(LUNOUT,*)'#RC+Born c/s per event: ',CS_EVT WRITE(LUNOUT,*)'#RC+Born events: ',NEVTOT_CORR ENDIF C. Smooth results C. Let's skip smoothing => we need to keep errors correctly. WRITE(LUNOUT,*)'--------- RESULTS SAVING ----------' C. Save all results C vm_born C. ----- Events file format ----- ver.1, One number per line C. (int) IVER :file format version (currently 1, reserved) C. (int) fNDPARS :size of parameters array C. (float) fDPARS(NPARAMS) : array of parameters (cuts and other DJANGOH settings) C. (int) NPTSX,NPTSY,NMODEX,NMODEY C. (float) XPTMIN,XPTMAX,YPTMIN,YPTMAX C. (int) NEVTS_BORN,IEVTS_BORN(NPTSX*NPTSY) C. (int) NEVTS_CORR,IEVTS_CORR(NPTSX*NPTSY) C. (int) fNEVTS_AddN (number of additional events tables) C. (int) fNEVTS_Add(NPTSX*NPTSY*fNEVTS_AddN) (additional data) C. (int) NCTRLSUM: currently NPTSX*NPTSY IVER=1 OPEN(NDATFILE,FILE='djangoh_r.evts',ERR=338,STATUS='REPLACE') WRITE(NDATFILE,*) IVER C. 2 additional parameters: FBINWX,FBINWY,ICUT,XMIN,XMAX,YMIN,YMAX,Q2MIN,Q2MAX,WMIN,S WRITE(NDATFILE,*) 13 WRITE(NDATFILE,491) FBINWX WRITE(NDATFILE,491) FBINWY WRITE(NDATFILE,*) ICUT WRITE(NDATFILE,491) XMIN WRITE(NDATFILE,491) XMAX WRITE(NDATFILE,491) YMIN WRITE(NDATFILE,491) YMAX WRITE(NDATFILE,491) Q2MIN WRITE(NDATFILE,491) Q2MAX WRITE(NDATFILE,491) WMIN WRITE(NDATFILE,492) GSP C. Save CS_EVT with 9-digits precision WRITE(NDATFILE,492) CS_EVT WRITE(NDATFILE,492) CS_EVT_BORN C. Point tables data WRITE(NDATFILE,*) NPTSX WRITE(NDATFILE,*) NPTSY WRITE(NDATFILE,*) NMODEX WRITE(NDATFILE,*) NMODEY WRITE(NDATFILE,491) XPTMIN WRITE(NDATFILE,491) XPTMAX WRITE(NDATFILE,491) YPTMIN WRITE(NDATFILE,491) YPTMAX WRITE(NDATFILE,*) NEVTOT_BORN DO 332 I1=1,NPTSX DO 332 I2=1,NPTSY WRITE(NDATFILE,*)NEVTS_BORN(I1,I2) 332 CONTINUE WRITE(NDATFILE,*)NEVTOT_CORR DO 333 I1=1,NPTSX DO 333 I2=1,NPTSY WRITE(NDATFILE,*)NEVTS_CORR(I1,I2) 333 CONTINUE C. Save additional event tables (if exist) WRITE(NDATFILE,*)NADDS IF(NADDS.GT.0) THEN DO 334 I3=1,NADDS DO 334 I1=1,NPTSX DO 334 I2=1,NPTSY WRITE(NDATFILE,*)NEVTS_ADD(I1,I2,I3) 334 CONTINUE ENDIF C. Save control sum (currently =NPTSX*NPTSY) NCTRLSUM=NPTSX*NPTSY WRITE(NDATFILE,*)NCTRLSUM CLOSE(NDATFILE) WRITE(LUNOUT,*)'djangoh_r.evts file saved' GOTO 339 338 WRITE(LUNOUT,*)'**** Error saving file! ****' 339 CONTINUE C. IF ((NEVTOT_CORR.GT.0).AND.(NEVTOT_BORN.GT.0)) THEN IF (NEVTOT_CORR.GT.0) THEN I1=IUPD_DATA_FILE() WRITE(LUNOUT,*)'+----------------------------------+' IF(I1.GE.3)THEN WRITE(LUNOUT,*)'| New data file created / replaced |' ELSEIF(I1.GE.2)THEN WRITE(LUNOUT,*)'| Data file updated |' ELSEIF(I1.GE.1)THEN WRITE(LUNOUT,*)'| Data file NOT updated |' WRITE(LUNOUT,*)'| No better precision achieved |' ELSE WRITE(LUNOUT,*)'| Error saving data file |' ENDIF WRITE(LUNOUT,*)'+----------------------------------+' ELSE WRITE(LUNOUT,*)'Results saving skipped' ENDIF C.********************************************************************** C. Save BORN cross section table for debug purposes * C. File format * C. #(int) IVER: file id/version (reserved, currently 1) * C. #(int) S: S-Mp2-Ml2 * C. #(int) NPTSX,NPTSY,NMODEX,NMODEY * C. #(float) XPTMIN,XPTMAX,YPTMIN,YPTMAX * C. #(float) DSIGMA/DXDY x(NPTSX*NPTSY) in [pb] * C. #(int) NCTRLSUM: currently NPTSX*NPTSY * C.********************************************************************** OPEN(NDATFILE,FILE='djangoh_r.born',ERR=378,STATUS='REPLACE') WRITE(NDATFILE,*) 1 WRITE(NDATFILE,492) GSP WRITE(NDATFILE,*) NPTSX WRITE(NDATFILE,*) NPTSY WRITE(NDATFILE,*) NMODEX WRITE(NDATFILE,*) NMODEY WRITE(NDATFILE,491) XPTMIN WRITE(NDATFILE,491) XPTMAX WRITE(NDATFILE,491) YPTMIN WRITE(NDATFILE,491) YPTMAX DO 370 I1=1,NPTSX DO 370 I2=1,NPTSY TEMP1=RPTSX(I1) TEMP2=RPTSY(I2) TEMP3=DJ_GET_BORN(TEMP1,TEMP2) IF(TEMP3.LE.0)TEMP3=-1D0 WRITE(NDATFILE,492)TEMP3 370 CONTINUE NCTRLSUM=NPTSX*NPTSY WRITE(NDATFILE,*)NCTRLSUM CLOSE(NDATFILE) WRITE(LUNOUT,*)'DEBUG: Born values file saved' GOTO 379 378 WRITE(LUNOUT,*)'DEBUG: Error saving Born value file!' 379 CONTINUE C. DEBUG variables log WRITE(LUNOUT,*)'NFBINFAILS=',NFBINFAILS(NADDS+1) DO 401 I3=1,NADDS WRITE(LUNOUT,*)'NFBINFAILS (ADD)=',NFBINFAILS(I3) 401 CONTINUE C...Original user procedure output IF (IHSONL.EQ.0) THEN WRITE(LUNOUT,3001) NTOT,NPASS,NSOPH,NFAILL,NFAILP IF (NTOT.NE.NPASS) WRITE(LUNOUT,3003) PARL(24) C...flavor distribution; new flavor coding (ckc) WRITE(LUNOUT,3004) IFLCNT(1),IFLCNT(2),IFLCNT(3) & ,IFLCNT(4),IFLCNT(5),IFLCNT(6) & ,IFLCNT(-1),IFLCNT(-2),IFLCNT(-3) & ,IFLCNT(-4),IFLCNT(-5),IFLCNT(-6) C...failures in hadronization WRITE(LUNOUT,3005) NFAILC(1),NFAILC(2),NFAILC(3) & ,NFAILC(4),NFAILC(5),NFAILC(6) & ,NFAILC(7),NFAILC(8),NFAILC(9) & ,NFAILC(10) ENDIF RETURN 491 FORMAT(E14.6) C 492 FORMAT(E17.9) 492 FORMAT(E20.12) 3001 FORMAT(/,' ******************************************************' F,'************************',/ F ,' Program performance ',/ F ,1X,I12,' Events were accepted by HERACLES,',/ F ,1X,I12,' Events passed fragmentation in LEPTO',/ F ,1X,I12,' Events passed fragmentation in SOPHIA,',/ F ,1X,I12,' Events failed fragmentation in LEPTO',/ F ,1X,I12,' Events failed fragmentation in SOPHIA') 3003 FORMAT(' Cross section was corrected: ',/ F ,' Not all events could be hadronized. ',/ F ,' Total cross section is now SIGTOT = ',E12.5,' pb',/) 3004 FORMAT(/,' Distribution of flavors (DJ):',/ F ,10X,'d',9X,'u',9X,'s',9X,'c',9X,'b',9X,'t',/,1X,6I10,/ F ,7X,'dbar',6X,'ubar',6X,'sbar',6X,'cbar',6X,'bbar',6X F ,'tbar',/,1X,6I10,/) 3005 FORMAT(/,' Errors in DJGEVT: Hadronization failed:',/ F ,' Nfail(1) = ',I8,' Nfail(2) = ',I8,/ F ,' Nfail(3) = ',I8,' Nfail(4) = ',I8,/ F ,' Nfail(5) = ',I8,' Nfail(6) = ',I8,/ F ,' Nfail(7) = ',I8,' Nfail(8) = ',I8,/ F ,' Nfail(9) = ',I8,' Nfail(10) = ',I8,/) END c*********************************************************************** c IDJR_FINDAFFBN: Find bins affected by {X,Y} event * c E.g. xleft=x/(1+width),... * c Note: every event may contribute several bins * c Uses bin width data: XPT_L,XPT_R,YPT_L,YPT_R * c*********************************************************************** FUNCTION IDJR_FINDAFFBN(X,Y,NBIN_XL,NBIN_XR,NBIN_YL,NBIN_YR) C IMPLICIT DOUBLE PRECISION (A-H,M,O-Z) IMPLICIT NONE DOUBLE PRECISION X,Y INTEGER IDJR_FINDAFFBN,NBIN_XL,NBIN_YL,NBIN_XR,NBIN_YR C. Copy of COMMON/DJRES_DATA1/ INTEGER NPTSX,NPTSY,NPTS_MAX,NPTS_MAX2,NMODEX,NMODEY,NCTRLSUM PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,CS_EVT,CS_EVT_BORN INTEGER NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN,NEVTS_CORR COMMON/DJRES_DATA1/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN(NPTS_MAX,NPTS_MAX) & ,NEVTS_CORR(NPTS_MAX,NPTS_MAX),CS_EVT,CS_EVT_BORN C. Copy of COMMON/DJRES_BINS/: bins info common block REAL*8 FBINWX,FBINWY,XPT_L,XPT_R,XBIN_W,YPT_L,YPT_R,YBIN_W COMMON/DJRES_BINS/FBINWX,FBINWY,XPT_L(NPTS_MAX),XPT_R(NPTS_MAX), & XBIN_W(NPTS_MAX),YPT_L(NPTS_MAX),YPT_R(NPTS_MAX),YBIN_W(NPTS_MAX) C. Locals INTEGER I1,I2 C... Find the most left/right bins for {X,Y} C... Returns 'bin-exist' flag (1 or 0) C... For faster bin search we presume: C... 1. The width of X- and Y- bins are independent. C... 2. For every bins n,n-1: Left(n-1) NXLEFT I1=1 I2=NPTSX 520 NXLEFT=I1+(I2-I1)/2 IF(RPTSX(NXLEFT).LT.X)THEN I1=NXLEFT IF(I2.GT.(I1+1))GOTO 520 C. Additional check for X==XPTMAX case if(RPTSX(I2).EQ.X)NXLEFT=I2 ELSEIF(RPTSX(NXLEFT).GT.X)THEN I2=NXLEFT IF(I2.GT.(I1+1))GOTO 520 NXLEFT=I1 ENDIF C. Find left Y-point in RPTSY set -> NYLEFT I1=1 I2=NPTSY 521 NYLEFT=I1+(I2-I1)/2 IF(RPTSY(NYLEFT).LT.Y)THEN I1=NYLEFT IF(I2.GT.(I1+1))GOTO 521 C. Additional check for Y==YPTMAX case if(RPTSY(I2).EQ.Y)NYLEFT=I2 ELSEIF(RPTSY(NYLEFT).GT.Y)THEN I2=NYLEFT IF(I2.GT.(I1+1))GOTO 521 NYLEFT=I1 ENDIF IDJR_FIND_LBIN=1 RETURN END c*********************************************************************** c DJ_GET_BORN: Returns Born cross section from DJANGOH (in pb) * c*********************************************************************** DOUBLE PRECISION FUNCTION DJ_GET_BORN(X,Y) IMPLICIT DOUBLE PRECISION (A-H,M,O-Z) DOUBLE PRECISION X,Y C. --------- DJANGOH declarations --------- C. declarations for heracles C ICHNN: event type, born(2) or brems(12) COMMON /HSCHNN/ ICHNN COMMON /HSPARL/ LPAR(20),LPARIN(12) COMMON /HSPARM/ POLARI,LLEPT,LQUA C. Locals DOUBLE PRECISION LPAR2_PREV LPAR2_PREV=LPAR(2) LPAR(2)=0 IF ((ICHNN.EQ.2).OR.(ICHNN.EQ.12)) THEN DJ_GET_BORN = HSSGCC(X,Y,LLEPT,POLARI,LQUA) ELSE C PRINT *,'ERROR: UNEXPECTED NC CALL' C STOP 8 DJ_GET_BORN = HSSGNC(X,Y,LLEPT,POLARI,LQUA) ENDIF LPAR(2)=LPAR2_PREV C. Convert [nb]->[pb] DJ_GET_BORN=DJ_GET_BORN*1D3 RETURN END c*********************************************************************** c IDJRES_PT_DELTA: Get relative correction value & calculate errors * c * c Input: NX,NY(point number) * c Result: DEL,ERR_STA,ERR_SYS * c. Returns: 1(ok), 0(no data) * c*********************************************************************** FUNCTION IDJRES_PT_DELTA(NX,NY,DEL,ERR_STA,ERR_SYS) IMPLICIT NONE INTEGER IDJRES_PT_DELTA,NX,NY,ITBL DOUBLE PRECISION DEL,ERR_STA,ERR_SYS C. Copy of COMMON/DJRES_DATA1/ INTEGER NPTSX,NPTSY,NPTS_MAX,NPTS_MAX2,NMODEX,NMODEY,NCTRLSUM PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,CS_EVT,CS_EVT_BORN INTEGER NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN,NEVTS_CORR COMMON/DJRES_DATA1/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN(NPTS_MAX,NPTS_MAX) & ,NEVTS_CORR(NPTS_MAX,NPTS_MAX),CS_EVT,CS_EVT_BORN C. Copy of COMMON/DJRES_BINS/: bins info common block REAL*8 FBINWX,FBINWY,XPT_L,XPT_R,XBIN_W,YPT_L,YPT_R,YBIN_W COMMON/DJRES_BINS/FBINWX,FBINWY,XPT_L(NPTS_MAX),XPT_R(NPTS_MAX), & XBIN_W(NPTS_MAX),YPT_L(NPTS_MAX),YPT_R(NPTS_MAX),YBIN_W(NPTS_MAX) C. declarations for heracles INTEGER LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND COMMON /HSUNTS/ LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND INTEGER INT2,INT3,ISAM2,ISAM3,IOPLOT,IPRINT,ICUT COMMON /HSOPTN/ INT2(5),INT3(15),ISAM2(5),ISAM3(15), * IOPLOT,IPRINT,ICUT INTEGER NEVENT,NEVE DOUBLE PRECISION SIGTOT,SIGTRR,SIGG,SIGGRR COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20) DOUBLE PRECISION SP,EELE,PELE,EPRO,PPRO COMMON /HSELAB/ SP,EELE,PELE,EPRO,PPRO DOUBLE PRECISION XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN COMMON /HSCUTS/ XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN DOUBLE PRECISION MEI,MEF,MQI,MQF,MEI2,MEF2,MQI2,MQF2,MPRO,MPRO2 COMMON /HSGSW1/ MEI,MEF,MQI,MQF,MEI2,MEF2,MQI2,MQF2,MPRO,MPRO2 C. Functions used INTEGER IDJR_FIND_LBIN DOUBLE PRECISION DJ_GET_BORN C. Locals INTEGER NCORR,NBORN,ITMP,ITMP2 INTEGER NPTXL,NPTXR,NPTYL,NPTYR,NEVXL,NEVXR,NEVYL,NEVYR REAL*8 TMP1 REAL*8 SIG_STA2,SIG_SYS2,ERR_LINX,ERR_LINY,ERR_LIN REAL*8 GSP,X,Y,Q2 REAL*8 DX,DY,X2,X1,Y2,Y1 C. Checks IDJRES_PT_DELTA=0 IF((NX.LT.1).OR.(NY.GT.NPTSX).OR.(NY.LT.1).OR.(NY.GT.NPTSY))RETURN IF(NEVTOT_CORR.LE.0)RETURN C. TODO: insert more checks for /DJRES_DATA1/ data here. GSP=SP-MPRO2-MEI2 X=RPTSX(NX) Y=RPTSY(NY) Q2=GSP*X*Y NCORR=NEVTS_CORR(NX,NY) IF(NCORR.LE.0)RETURN IF((X.LE.0D0).OR.(Y.LE.0D0))RETURN C. Get delta C. Old code for delta via separate RC and BORN samplings C. NBORN=NEVTS_BORN(NX,NY) C. IF((NBORN.LE.0).OR.(CS_EVT_BORN.LE.0).OR.(CS_EVT.LE.0))RETURN C. DEL=NCORR C. IF(NBORN.GT.0)DEL=DEL*CS_EVT/DFLOAT(NBORN)/CS_EVT_BORN-1D0 C. New code for delta via RC sampling and internal BORN routine IF((FBINWX.LE.0).OR.(FBINWY.LE.0))RETURN X2=X*(1D0+FBINWX) X1=X*(1D0-FBINWX) Y2=Y*(1D0+FBINWY) Y1=Y*(1D0-FBINWY) C. Check for borders vicinity IF(X2.GT.XMAX)X2=XMAX IF(X1.LT.XMIN)X1=XMIN C.Note: We don't need to check WMIN cut (ICUT.GT.2) IF(ICUT.EQ.3)THEN IF(Y1.LT.YMIN)Y1=YMIN C.Check for YMAX is performed in IUPD_DATA_FILE ==> next check is superfluous IF(Y2.GT.YMAX)Y2=YMAX IF((GSP*X2*Y2).GT.Q2MAX)RETURN ENDIF C.Check for Q2min is performed in IUPD_DATA_FILE ==> next check is superfluous IF((GSP*X1*Y1).LT.Q2MIN)RETURN DX=X2-X1 DY=Y2-Y1 IF((DX.LE.0D0).OR.(DY.LE.0D0))RETURN TMP1=DJ_GET_BORN(X,Y) IF(TMP1.EQ.0D0)RETURN DEL=DFLOAT(NCORR)*CS_EVT/DX/DY/TMP1-1D0 C. Calculate error SIG_STA2=0D0 SIG_SYS2=0D0 C. Statistical error C. Nt: SigmaC=Corr*Sqrt(1/Nbin-1/Nevts)~Corr*Sqrt(1/Nbin) C. Nt: SigmaB=Born*Sqrt(1/Nbin-1/Nevts)~Born*Sqrt(1/Nbin) SIG_STA2=SIG_STA2+1D0/DFLOAT(NCORR)-1D0/DFLOAT(NEVTOT_CORR) C. Systematical error (due to finite bin width) NEVXL=0 NEVXR=0 NEVYL=0 NEVYR=0 ITMP=IDJR_FIND_LBIN(X*(1D0-FBINWX),Y,NPTXL,ITMP2) IF(ITMP.GT.0)THEN NEVXL=NEVTS_CORR(NPTXL,NY) ENDIF ITMP=IDJR_FIND_LBIN(X,Y*(1D0-FBINWY),ITMP2,NPTYL) IF(ITMP.GT.0)THEN NEVYL=NEVTS_CORR(NX,NPTYL) ENDIF ITMP=IDJR_FIND_LBIN(X*(1D0+FBINWX),Y,NPTXR,ITMP2) IF(ITMP.GT.0)THEN NEVXR=NEVTS_CORR(NPTXR,NY) ENDIF ITMP=IDJR_FIND_LBIN(X,Y*(1D0+FBINWY),ITMP2,NPTYR) IF(ITMP.GT.0)THEN NEVYR=NEVTS_CORR(NX,NPTYR) ENDIF ERR_LINX=0D0 IF((NEVXL.GT.0).AND.(NEVXR.GT.0))THEN ERR_LINX=DFLOAT(NEVXL+NEVXR-2*NCORR)/2D0 ELSEIF(NEVXL.GT.0)THEN ERR_LINX=DFLOAT(NEVXL-NCORR) ELSEIF(NEVXR.GT.0)THEN ERR_LINX=DFLOAT(NEVXR-NCORR) ELSE ERR_LINX=0D0 ENDIF ERR_LINY=0D0 IF((NEVYL.GT.0).AND.(NEVYR.GT.0))THEN ERR_LINY=DFLOAT(NEVYL+NEVYR-2*NCORR)/2D0 ELSEIF(NEVYL.GT.0)THEN ERR_LINY=DFLOAT(NEVYL-NCORR) ELSEIF(NEVYR.GT.0)THEN ERR_LINY=DFLOAT(NEVYR-NCORR) ELSE ERR_LINY=0D0 ENDIF C. IF(ERR_LINX.LT.0)ERR_LINX=-ERR_LINX C. IF(ERR_LINY.LT.0)ERR_LINY=-ERR_LINY C ERR_LIN=(ERR_LINX+ERR_LINY)/2/NCORR ERR_LIN=DMAX1(DABS(ERR_LINX),DABS(ERR_LINY)) ERR_LIN=ERR_LIN/NCORR IF(ERR_LIN.LT.0D0)ERR_LIN=-ERR_LIN IF(ERR_LIN.GT.1D0)ERR_LIN=1D0 SIG_SYS2=SIG_SYS2+ERR_LIN**2 C. Note: Born is presumed to be exact here. C. Total error C. Nt: Sigma{Delta}=Corr/Born*Sqrt((SigmaCorr/Corr)^2+(SigmaBorn/Born)^2) C. ERR_STA=DSQRT(SIG_STA2+...)*(DEL+1D0) ERR_STA=DSQRT(SIG_STA2)*(DEL+1D0) C. ERR_SYS=DSQRT(SIG_SYS2+...)*(DEL+1D0) ERR_SYS=DSQRT(SIG_SYS2)*(DEL+1D0) IDJRES_PT_DELTA=1 RETURN END c*********************************************************************** c IUPD_DATA_FILE: Update results file * c * c. File format * c. #(int) IVER: file version (reserved, currently 2) * c. #(int) I2ERRS: 1(0) the stat. & sys. errors are separated * c. #(int) S: S-Mp2-Ml2 * c. #(int) NPTSX,NPTSY,NMODEX,NMODEY * c. #(float) XPTMIN,XPTMAX,YPTMIN,YPTMAX * c. #(int) DELTA,ERRSTA[,ERRSYS] x(NPTSX*NPTSY) * c. #(int) NCTRLSUM: currently NPTSX*NPTSY * c*********************************************************************** FUNCTION IUPD_DATA_FILE() IMPLICIT NONE INTEGER IUPD_DATA_FILE C. Copy of COMMON/DJRES_DATA1/ INTEGER NPTSX,NPTSY,NPTS_MAX,NPTS_MAX2,NMODEX,NMODEY,NCTRLSUM PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,CS_EVT,CS_EVT_BORN INTEGER NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN,NEVTS_CORR COMMON/DJRES_DATA1/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS,NEVTOT_BORN,NEVTOT_CORR,NEVTS_BORN(NPTS_MAX,NPTS_MAX) & ,NEVTS_CORR(NPTS_MAX,NPTS_MAX),CS_EVT,CS_EVT_BORN C. Copy of COMMON/DJRES_BINS/: bins info common block REAL*8 FBINWX,FBINWY,XPT_L,XPT_R,XBIN_W,YPT_L,YPT_R,YBIN_W COMMON/DJRES_BINS/FBINWX,FBINWY,XPT_L(NPTS_MAX),XPT_R(NPTS_MAX), & XBIN_W(NPTS_MAX),YPT_L(NPTS_MAX),YPT_R(NPTS_MAX),YBIN_W(NPTS_MAX) C. Results array REAL*8 DELTA(NPTS_MAX,NPTS_MAX) REAL*8 ERRSTA(NPTS_MAX,NPTS_MAX),ERRSYS(NPTS_MAX,NPTS_MAX) REAL*8 VAL_NO_DATA,ERR_NO_DATA PARAMETER(VAL_NO_DATA=1D6, ERR_NO_DATA=1D0) C. declarations for heracles INTEGER LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND COMMON /HSUNTS/ LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND INTEGER INT2,INT3,ISAM2,ISAM3,IOPLOT,IPRINT,ICUT COMMON /HSOPTN/ INT2(5),INT3(15),ISAM2(5),ISAM3(15), * IOPLOT,IPRINT,ICUT INTEGER NEVENT,NEVE DOUBLE PRECISION SIGTOT,SIGTRR,SIGG,SIGGRR COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20) DOUBLE PRECISION SP,EELE,PELE,EPRO,PPRO COMMON /HSELAB/ SP,EELE,PELE,EPRO,PPRO DOUBLE PRECISION XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN COMMON /HSCUTS/ XMIN,XMAX,Q2MIN,Q2MAX,YMIN,YMAX,WMIN,GMIN DOUBLE PRECISION MEI,MEF,MQI,MQF,MEI2,MEF2,MQI2,MQF2,MPRO,MPRO2 COMMON /HSGSW1/ MEI,MEF,MQI,MQF,MEI2,MEF2,MQI2,MQF2,MPRO,MPRO2 C. Functions used INTEGER IDJRES_PT_DELTA C. Locals INTEGER NRESFILE PARAMETER (NRESFILE=32) INTEGER I1,I2,ITMP1,ITMP2,ITMP3,ITMP4,IRESLOADED,NNEWVALS REAL*8 TMP1,TMP2,TMP3,TMP4 REAL*8 GSP,X,Y,Q2 REAL*8 DEL_NEW,ERR_NEW,ERR_NEW1,ERR_NEW2,ERR_OLD INTEGER ICUTS_OK INTEGER IVER,I2ERRS C. IUPD_DATA_FILE=0 IRESLOADED=0 I2ERRS=2 TMP1=0 TMP2=0 TMP3=0 TMP4=0 IVER=0 GSP=SP-MPRO2-MEI2 C. Try to load previous file OPEN(NRESFILE,FILE='djangoh_r.dat',ERR=418,STATUS='OLD') WRITE(LUNOUT,*)'Loading previous results file' READ(NRESFILE,*)IVER IF (IVER.EQ.2) THEN IRESLOADED=1 READ(NRESFILE,*)I2ERRS IF((I2ERRS.GT.1).OR.(I2ERRS.LT.0))THEN WRITE(LUNOUT,*)'Unexpected I2ERRS value: ',I2ERRS IRESLOADED=0 ENDIF READ(NRESFILE,*)TMP1 C. Compare these parameters to current ones IF(DABS(TMP1-GSP).GT.(GSP*1D-5))THEN WRITE(LUNOUT,*)'Different energy in old file: ',TMP1,GSP IRESLOADED=0 ENDIF IF(IRESLOADED.EQ.1)THEN READ(NRESFILE,*) ITMP1 READ(NRESFILE,*) ITMP2 READ(NRESFILE,*) ITMP3 READ(NRESFILE,*) ITMP4 IF ((ITMP1.NE.NPTSX).OR.(ITMP2.NE.NPTSY).OR. & (ITMP3.NE.NMODEX).OR.(ITMP4.NE.NMODEY)) THEN WRITE(LUNOUT,*)'X,Y grid parameters are different' WRITE(LUNOUT,*)'NPTSX (old,new):',ITMP1,NPTSX WRITE(LUNOUT,*)'NPTSY (old,new):',ITMP2,NPTSY WRITE(LUNOUT,*)'NMODEX (old,new):',ITMP3,NMODEX WRITE(LUNOUT,*)'NMODEY (old,new):',ITMP4,NMODEY IRESLOADED=0 ENDIF ENDIF IF(IRESLOADED.EQ.1)THEN READ(NRESFILE,*) TMP1 READ(NRESFILE,*) TMP2 READ(NRESFILE,*) TMP3 READ(NRESFILE,*) TMP4 IF ((TMP1.NE.XPTMIN).OR.(TMP2.NE.XPTMAX).OR. & (TMP3.NE.YPTMIN).OR.(TMP4.NE.YPTMAX)) THEN WRITE(LUNOUT,*)'X,Y margins are different' WRITE(LUNOUT,*)'XPTMIN (old,new):',TMP1,XPTMIN WRITE(LUNOUT,*)'XPTMAX (old,new):',TMP2,XPTMAX WRITE(LUNOUT,*)'YPTMIN (old,new):',TMP3,YPTMIN WRITE(LUNOUT,*)'YPTMAX (old,new):',TMP4,YPTMAX IRESLOADED=0 ENDIF ENDIF IF (IRESLOADED.EQ.1) THEN DO 412 I1=1,NPTSX DO 412 I2=1,NPTSY READ(NRESFILE,*)TMP1 DELTA(I1,I2)=TMP1 READ(NRESFILE,*)TMP2 ERRSTA(I1,I2)=TMP2 IF(I2ERRS.EQ.1)THEN READ(NRESFILE,*)TMP3 ERRSYS(I1,I2)=TMP3 ELSE ERRSYS(I1,I2)=0D0 ENDIF 412 CONTINUE READ(NRESFILE,*)NCTRLSUM IF (NCTRLSUM.NE.(NPTSX*NPTSY)) THEN WRITE(LUNOUT,*)'Invalid control sum read:',NCTRLSUM IRESLOADED=0 ENDIF ENDIF ELSE WRITE(LUNOUT,*)'Unexpected version of DAT-file' WRITE(LUNOUT,*)'The previous results will be overwritten' IRESLOADED=0 ENDIF CLOSE(NRESFILE) IF(IRESLOADED.EQ.1)THEN WRITE(LUNOUT,*)'Previous results file succesfully loaded' ELSE WRITE(LUNOUT,*)'Results file loading failed' ENDIF GOTO 419 418 WRITE(LUNOUT,*)'No results file exist yet' IRESLOADED=0 419 CONTINUE C. Prepate 'no-data' values if no data read IF (IRESLOADED.EQ.0) THEN DO 421 I1=1,NPTSX DO 421 I2=1,NPTSY DELTA(I1,I2)=VAL_NO_DATA ERRSTA(I1,I2)=ERR_NO_DATA ERRSYS(I1,I2)=ERR_NO_DATA 421 CONTINUE I2ERRS=1 IVER=2 ENDIF C. Update results WRITE(LUNOUT,*)'Updating results...' NNEWVALS=0 DO 431 I1=1,NPTSX DO 431 I2=1,NPTSY DEL_NEW=VAL_NO_DATA ERR_NEW=ERR_NO_DATA ERR_NEW1=ERR_NO_DATA ERR_NEW2=ERR_NO_DATA ICUTS_OK=1 X=RPTSX(I1) Y=RPTSY(I2) Q2=GSP*X*Y IF((Q2*(1D0-FBINWX)*(1D0-FBINWY)).LT.Q2MIN)ICUTS_OK=0; IF((X*(1D0+FBINWX)).GT.XMAX)ICUTS_OK=0; IF(ICUT.EQ.3)THEN IF((Y*(1D0-FBINWY)).LT.YMIN)ICUTS_OK=0; ENDIF C. Note: The complete check of cuts is implemented in IDJRES_PT_DELTA() IF(ICUTS_OK.GT.0)THEN ITMP1=IDJRES_PT_DELTA(I1,I2,DEL_NEW,ERR_NEW1,ERR_NEW2) IF(ITMP1.GT.0)THEN ERR_NEW=DSQRT(ERR_NEW1**2+ERR_NEW2**2) ELSE DEL_NEW=VAL_NO_DATA ERR_NEW=ERR_NO_DATA ERR_NEW1=ERR_NO_DATA ERR_NEW2=ERR_NO_DATA ENDIF ENDIF ERR_OLD=ERR_NO_DATA IF ((ERRSYS(I1,I2).LT.ERR_NO_DATA).AND & .(ERRSTA(I1,I2).LT.ERR_NO_DATA)) THEN IF(I2ERRS.EQ.1)THEN ERR_OLD=DSQRT(ERRSTA(I1,I2)**2+ERRSYS(I1,I2)**2) ELSE ERR_OLD=ERRSTA(I1,I2) ENDIF ENDIF IF(ERR_NEW.LT.ERR_OLD)THEN DELTA(I1,I2)=DEL_NEW ERRSTA(I1,I2)=ERR_NEW1 ERRSYS(I1,I2)=ERR_NEW2 NNEWVALS=NNEWVALS+1 ENDIF 431 CONTINUE WRITE(LUNOUT,*)'Results updated. Updated values: ',NNEWVALS C. Return if no new results IF((NNEWVALS.LE.0).AND.(IRESLOADED.GT.0))THEN IUPD_DATA_FILE=1 RETURN ENDIF C. Save results OPEN(NRESFILE,FILE='djangoh_r.dat',ERR=478,STATUS='REPLACE') WRITE(NRESFILE,*) IVER WRITE(NRESFILE,*) I2ERRS WRITE(NRESFILE,492) GSP C. Point tables data WRITE(NRESFILE,*) NPTSX WRITE(NRESFILE,*) NPTSY WRITE(NRESFILE,*) NMODEX WRITE(NRESFILE,*) NMODEY WRITE(NRESFILE,491) XPTMIN WRITE(NRESFILE,491) XPTMAX WRITE(NRESFILE,491) YPTMIN WRITE(NRESFILE,491) YPTMAX DO 472 I1=1,NPTSX DO 472 I2=1,NPTSY WRITE(NRESFILE,491)DELTA(I1,I2) WRITE(NRESFILE,491)ERRSTA(I1,I2) IF(I2ERRS.EQ.1)WRITE(NRESFILE,491)ERRSYS(I1,I2) 472 CONTINUE C. Save control sum (currently =NPTSX*NPTSY) NCTRLSUM=NPTSX*NPTSY WRITE(NRESFILE,*)NCTRLSUM CLOSE(NRESFILE) WRITE(LUNOUT,*)'djangoh_r.dat file saved' IF(NNEWVALS.LE.0)THEN WRITE(LUNOUT,*)'No valid results exist. Empty file is saved.' IUPD_DATA_FILE=1 ELSE IUPD_DATA_FILE=3 IF (IRESLOADED.EQ.1)IUPD_DATA_FILE=2 ENDIF GOTO 479 478 WRITE(LUNOUT,*)'**** Error saving file! ****' IUPD_DATA_FILE=0 479 CONTINUE RETURN 491 FORMAT(E14.6) C 492 FORMAT(E17.9) 492 FORMAT(E20.12) END