* Example for RUN SUBROUTINE RUNTST * * This is the standard RUN test routine: call it in MAIN * needs scratch files 51 and 52 * generates DATA and MOCA and BACK ntuple data * call the unfolding * ------------------------------------------------------- CHARACTER*4 NAME C CHARACTER*80 STRING C CALL GETARG(IARGC(),STRING) C WRITE(*,*) STRING * ... CALL CAIM('DATA DATA 5000 ') CALL CAIM('MONTECARLO MOCA 20000.0 ') CALL CAIM('BACKGROUND BACK 500.0 ') C CALL CAIM('DATA DATA ') C CALL CAIM('MONTECARLO MOCA ') C CALL CAIM('BACKGROUND BACK 500.0 ') CALL CAIM('VARIABLE 0 ! Unfolding variable X') C CALL CAIM('VARIABLE 1 120 ! Measured variable Y(1)') C CALL CAIM('VARIABLE 2 10 ! Measured variable Y(2)') CALL CAIM('VARIABLE 1 40 ! Measured variable Y(1)') CALL CAIM('VARIABLE 2 ! Measured variable Y(2)') CALL CAIM('TITLE ! Test of the unfolding program') CALL CAIM('KNOTS 25') CALL CAIM('SMOOTHX') CALL CAIM('NRDF 10') CALL CAIM('XLIMITS 0.0 2.0') CALL CAIM('LINWEIGHT 0.1') CALL CAIM('XBINS 0.0 0.16 0.32 0.48 0.66 0.84 1.01 1.19 + 1.36 1.54 1.71 1.87 2.0') CALL CAIM('PRINTFLAG 2') C CALL CAIM('CONSTRAINT 1.2 0.25') C CALL CAIM('CONSTRAINT 1.4 0.95') C CALL CAIM('KNOTS 43') C CALL CAIM('NRDF 20') C CALL CAIM('XLIMITS 0.0 2.0') C CALL CAIM('LINWEIGHT 0.1') * generate ntuples ... CALL RUTEST(RATLUM) * ... getting back ratio of data:moca luminosity 10 CALL RPASS1(NAME) IF(NAME.EQ.' ') GOTO 20 GOTO 10 20 CONTINUE CALL RPASS2 30 CALL RPASS3(NAME) * check end of data IF(NAME.EQ.' ') GOTO 40 * ... GOTO 30 40 CONTINUE * ... and printer plots CALL RUNPRT 100 RETURN END SUBROUTINE RUTEST(RATLUM) * Standard test subroutine for RUN * * generate test sample of data and moca n-tuples * * common for TEST histograms --------------------------------------- PARAMETER (NBINS=80,NDIM=10+NBINS) COMMON/CRTEST/FTI(101), + DATAXT(NDIM),DATAXA(NDIM),DATAYT(NDIM),DATAYA(NDIM),DATAYB(NDIM), + RANDXT(NDIM),RANDXA(NDIM),RANDYT(NDIM),RANDYA(NDIM),BACKYT(NDIM), + NDATA,NMOCA,NDACK,NBACK SAVE /CRTEST/ * end of common for TEST histograms -------------------------------- REAL X,Y(4) REAL DLUMI,MLUMI,RATLUM * Parameters and functions for Test ntuple generation -------------- * Parameters for FTX statement function DATA B1/1.0/,X1/0.4/,G1/2.0/, + B2/10.0/,X2/0.8/,G2/0.2/, + B3/ 5.0/,X3/1.5/,G3/0.2/ * Resolution DATA SIGMA/0.1/ * Statement functions FTX(X)=B1*G1**2/((X-X1)**2+G1**2) + +B2*G2**2/((X-X2)**2+G2**2) + +B3*G3**2/((X-X3)**2+G3**2) * Acceptance function ACC(X)=1.0-0.50*(X-1.00)**2 * Transformation function TRF(X)=X-0.050*X**2 * ... * number of events to be generated NDATA= 5000 ! number of data events NMOCA=20000 ! number of moca ntuples NDACK= 500 ! number of backgrund events within data NBACK= 5000 ! background ntuples * end of: Parameters and functions for Test ntuple generation ------ * Version for LOOK ntuples - booking * CALL BVEC( 77, 1,4) * CALL BVEC( 77, 2,4) * CALL BVEC( 77, 3,4) * initialize histograms for true and measured distribution CALL AHISTS(DATAXT,NBINS,0.0,2.0) ! all data X CALL AHISTS(DATAXA,NBINS,0.0,2.0) ! accepted data X CALL AHISTS(DATAYT,NBINS,0.0,2.0) ! all data Y CALL AHISTS(DATAYA,NBINS,0.0,2.0) ! accepted data Y CALL AHISTS(DATAYB,NBINS,0.0,2.0) ! backgrnd data Y CALL AHISTS(RANDXT,NBINS,0.0,2.0) ! all moca X CALL AHISTS(RANDXA,NBINS,0.0,2.0) ! accepted moca X CALL AHISTS(RANDYT,NBINS,0.0,2.0) ! all moca Y CALL AHISTS(RANDYA,NBINS,0.0,2.0) ! accepted moca Y CALL AHISTS(BACKYT,NBINS,0.0,2.0) ! all back Y CALL SMTEXT('RUTEST: Standard RUN test subroutine') * ------------------------------------------ * find maximum of true function and integral XA=0.0 XB=2.0 FINT=0.0 DO L=1,5 FMAX=0.0 DO I=1,101 * x value X=0.01*(XA*FLOAT(101-I)+XB*FLOAT(I-1)) IF(L.EQ.1) THEN * Integration (for first loop L=1) FAC=0.02 IF(I.EQ.1.OR.I.EQ.101) FAC=0.01 FTXX=FTX(X) FINT=FINT+FTXX*FAC * ... for function plot FTI(I)=FTXX END IF IF(FTX(X).GT.FMAX) THEN * maximum of true function XMAX=X FMAX=FTX(X) END IF END DO C WRITE(*,*) L,XA,XB C WRITE(6,101) FINT,XMAX * reduce range of maximum by factor 0.02 DX=XB-XA XA=XMAX-0.01*DX XB=XMAX+0.01*DX END DO * normalize function to unit integral DO I=1,101 FTI(I)=FTI(I)/FINT END DO WRITE(6,101) FINT,XMAX CALL BGTEXT('GENERATION OF ') CALL BGTEXT(' TEST RECORDS') DO IDM=1,3 * ntuples * IDM = 1 Data * IDM = 2 Moca * IDM = 3 background IF(IDM.EQ.1) NN=NDATA IF(IDM.EQ.2) NN=NMOCA IF(IDM.EQ.3) NN=NDACK+NBACK * generate NN ntuples DO N=1,NN * weight initially allways 1 WT=1.0 * generate x, uniformly distributed between 0 and 2 30 RANDOM=RNUM() C IF(N.LE.50) WRITE(*,*) 'random ',N,RANDOM X=2.0*RANDOM * different treatment for different cases IF(IDM.EQ.1) THEN * data: hit-miss monte carlo to get distribution FTX IF(FMAX*RNUM().GT.FTX(X)) GOTO 30 * apply acceptance probability (random) - data lost IF(RNUM().GT.ACC(X)) WT=0.0 * apply systematic transformation and smearing - no loss 31 Y(1)=TRF(X)+SIGMA*RNGA() IF(Y(1).LT.0.0.OR.Y(1).GT.2.0) GOTO 31 * entries to histogram CALL AHISTU(DATAXT,X) ! all data X IF(WT.NE.0.0) CALL AHISTU(DATAXA,X) ! accepted data X CALL AHISTU(DATAYT,Y(1)) ! all data Y IF(WT.NE.0.0) CALL AHISTU(DATAYA,Y(1)) ! accepted data Y ELSE IF(IDM.EQ.2) THEN * Moca: uniform * apply acceptance probability (random) - ntuple unmeasurable RAC=RNUM() ACF=ACC(X) IF(RAC.GT.ACF) THEN WT=-ABS(WT) END IF 32 Y(1)=TRF(X)+SIGMA*RNGA() IF(Y(1).LT.0.0.OR.Y(1).GT.2.0) GOTO 32 * entries to weighted histograms CALL AHISTW(RANDXT,X,ABS(WT)) ! all moca X IF(WT.GT.0.0) CALL AHISTW(RANDXA,X,WT) ! accepted moca X CALL AHISTW(RANDYT,Y(1),ABS(WT)) ! all moca Y IF(WT.GT.0.0) CALL AHISTW(RANDYA,Y(1),WT) ! accepted moca Y ELSE IF(IDM.EQ.3) THEN * background: continuous distibution X=0.5*X*X Y(1)=X * apply acceptance probability for background too, but no loss IF(RNUM().GT.ACC(X)) GOTO 30 * no transformation and smearing for background Y(1)=X IF(N.LE.NDACK) THEN * data entries to histogram CALL AHISTU(DATAYT,Y(1)) ! all data Y CALL AHISTU(DATAYA,Y(1)) ! accepted data Y CALL AHISTU(DATAYB,Y(1)) ! backgrnd data Y ELSE * back entries to histogram CALL AHISTU(BACKYT,Y(1)) ! all back Y END IF END IF * second measured value with weak dependence on true X Y(2)=0.5*(1.0-SQRT(RNUM()))*X * store weight and true X in elements 3 and 4 Y(3)=WT Y(4)=X * write to ntuple file A and add to histogram IF(IDM.EQ.1.AND.WT.NE.0.0) THEN C CALL SVEC(77,1,Y) CALL UEVENT('DATA',2,1.0,0.0,Y) ELSE IF(IDM.EQ.2) THEN C CALL SVEC(77,2,Y) CALL UEVENT('MOCA',2,WT,X,Y) ELSE IF(IDM.EQ.3) THEN IF(N.LE.NDACK) THEN * add to data C CALL SVEC(77,1,Y) CALL UEVENT('DATA',2,1.0,0.0,Y) ELSE * add to background C CALL SVEC(77,3,Y) CALL UEVENT('BACK',2,1.0,0.0,Y) END IF END IF END DO END DO * close ntuple file A CALL UEVENT('DATA',-1,X,WT,Y) * integrated luminosities for data ... DLUMI=FLOAT(NDATA) ! assuming total cross section of 1 * ... and for Monte Carlo ... MLUMI=FLOAT(NMOCA) ! for constant function * ratio of data to monte Carlo luminosity is: RATLUM=DLUMI/MLUMI * printout CALL SMTEXT('Histogram of true events') CALL AHISTP(DATAXT) C WRITE(*,*) 'DATAXT ',DATAXT CALL SMTEXT('Histogram of measured events') CALL AHISTP(DATAYT) 100 RETURN 101 FORMAT(' Integral of true function is',G15.5, 1 ' and maximum is at',G15.5) END