subroutine chi(ntyp,nhist,iteration,nbinmax,nbins, |asiz,useOverflow,useMCerror,filename,mainid,outdir, | obgnorms,ncnorms,cohnorms,ccnorms,obg_count,istep) implicit none ! Input/Output: integer ntyp integer nhist integer iteration integer nbinmax integer nbins(3) integer asiz integer useOverflow integer useMCerror character*80 filename(15,3) integer mainid(22,3) character*45 outdir(22) double precision obgnorms(22,25) double precision ncnorms(22,0:25,0:4) double precision cohnorms(22,0:25,0:4) double precision ccnorms(22,0:25,0:4) double precision obg_count(22,15) double precision istep ! Paw Declarations: integer nhbook parameter (nhbook=300000000) real paw common / pawc / paw(nhbook) integer istat,icycle real HI,HIE ! Do loop variables integer ii,jj,kk integer ihist,ityp integer ifloat integer inorm ! Contents of histogram bins: double precision histvec(ntyp,nbinmax,3) double precision histerr(ntyp,nbinmax,3) ! Temporary vectors for extracting histograms ! (to be put into "histvec" and "histerr") real tempvecNC(nbins(1)) real tempvecCC(nbins(2)) real tempvecCoh(nbins(3)) ! Chisq. array double precision chisq(asiz) real chisqr(asiz) !--> Input to Hbook call integer hid(ntyp,nhist,3) !--> Histogram IDs integer nzerobins(3) !--> Count of bins with no data events with ! ifloat index double precision CohNorm,NCnorm,CCnorm double precision MaxNorm,MinNorm double precision norm double precision tnorm(2:5) !--> Temporary normalizations double precision minchi double precision chisqtmp double precision errortmp double precision sigma character*80 hfile integer length call HLIMIT(nhbook) !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@! ! Add constant to mainid for each component type so ! that histogram ID's do not overlap in memory. do ifloat=1,3 do ihist=1,nhist do ityp=1,ntyp hid(ityp,ihist,ifloat) = 100000000*ifloat+ | 10000000*ityp + | mainid(ihist,ifloat) enddo enddo enddo !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@! !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Open input hbook file and get histograms do ifloat=1,3 do ityp=1,ntyp call hropen(1,'input',filename(ityp,ifloat),' ',1024,istat) call hcdir('//input',' ') do ihist=1,nhist call hrin(mainid(ihist,ifloat),9999, | ifloat*100000000+ityp*10000000) enddo call hrend('input') close(1) enddo enddo !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ !|||||||||||||||||||||||||||||||||||||||||||||||||||! !+++++++++++++++++++++++++++++++++++++++++++++++++++! ! ! do ihist = 1, nhist !--> Loop over histogram sets ! ! ! !+++++++++++++++++++++++++++++++++++++++++++++++++++! !|||||||||||||||||||||||||||||||||||||||||||||||||||! !------------------------------------------------ ! Dump histogram contents into arrays do ityp=1,ntyp call hunpak(hid(ityp,ihist,1),tempvecNC,' ',0) call hunpak(hid(ityp,ihist,2),tempvecCC,' ',0) call hunpak(hid(ityp,ihist,3),tempvecCoh,' ',0) do ii=1,nbins(1) histvec(ityp,ii,1) = dble(tempvecNC(ii)) enddo do ii=1,nbins(2) histvec(ityp,ii,2) = dble(tempvecCC(ii)) enddo do ii=1,nbins(3) histvec(ityp,ii,3) = dble(tempvecCoh(ii)) enddo enddo !------------------------------------------------ ! Dump histogram errors into arrays ! Data do ifloat=1,3 do ii=1,nbins(ifloat) histerr(1,ii,ifloat) = dsqrt(histvec(1,ii,ifloat)) enddo enddo ! MC do ityp=2,ntyp call hunpke(hid(ityp,ihist,1),tempvecNC,' ',0) call hunpke(hid(ityp,ihist,2),tempvecCC,' ',0) call hunpke(hid(ityp,ihist,3),tempvecCoh,' ',0) do ii=1,nbins(1) histerr(ityp,ii,1) = dble(tempvecNC(ii)) enddo do ii=1,nbins(2) histerr(ityp,ii,2) = dble(tempvecCC(ii)) enddo do ii=1,nbins(3) histerr(ityp,ii,3) = dble(tempvecCoh(ii)) enddo enddo !------------------------------------------------ ! Adding overflow to last bin if desired: if (useOverflow.eq.1) then do ifloat=1,3 do ityp=1,ntyp histvec(ityp,nbins(ifloat),ifloat) = | histvec(ityp,nbins(ifloat),ifloat) + | dble(HI(hid(ityp,ihist,ifloat),nbins(ifloat)+1)) histerr(ityp,nbins(ifloat),ifloat) = | dsqrt(histerr(ityp,nbins(ifloat),ifloat)**2 + | dble(HIE(hid(ityp,ihist,ifloat),nbins(ifloat)+1))**2) enddo enddo endif !------------------------------------------------ !=========================================! ! Data Bin Value Checks ! length=LEN_TRIM(outdir(ihist)) ! hfile = 'output/'//outdir(ihist)(1:length)! | //'/warning.txt' ! open(55,file=hfile,status='UNKNOWN') ! do ifloat=1,3 ! nzerobins(ifloat) = 0 ! do ii=1,nbins(ifloat) ! if(histvec(1,ii,ifloat).eq.0.0D0) then ! write(55,*)'' ! write(55,*)'!!! Warning !!!' ! write(55,*)'Data in Bin is Zero' ! write(55,*)'Not using Bin Number: ',ii ! write(55,*)'Variable Number: ',ihist ! write(55,*)'Float type: ',ifloat ! write(55,*)'HistID: ', ! | mainid(ihist,ifloat) ! write(55,*)'' ! nzerobins(ifloat)= ! | nzerobins(ifloat) + 1 ! endif ! enddo ! enddo ! close(55) ! !=========================================! do ifloat=1,3 !--> Loop over floated MC compnents. ! 1 = NCDIS ! 2 = CCDIS ! 3 = CohRho0 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! Set current Norm Factors: if(ifloat.eq.1) then NCnorm = ncnorms(ihist,iteration-1,0) CCnorm = ccnorms(ihist,iteration-1,0) elseif(ifloat.eq.2) then NCnorm = ncnorms(ihist,iteration ,0) CCnorm = ccnorms(ihist,iteration-1,0) elseif(ifloat.eq.3) then NCnorm = ncnorms(ihist,iteration,0) CCnorm = ccnorms(ihist,iteration,0) endif CohNorm=cohnorms(ihist,iteration-1,0) !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx !00000000000000000000000000000000000000000000000000000000! ! Set OBG normalization factor: ! First remove floated components from data: obgnorms(ihist,iteration) = obg_count(ihist,1) - | CohNorm*obg_count(ihist,2) - | NCnorm*obg_count(ihist,3) - | CCnorm*obg_count(ihist,4) ! Remove remaining components: if(nhist.ge.6) then do ityp=6,ntyp obgnorms(ihist,iteration) = | obgnorms(ihist,iteration) - obg_count(ihist,ityp) enddo endif ! Now find OBG normalization: obgnorms(ihist,iteration) = | obgnorms(ihist,iteration)/obg_count(ihist,5) !00000000000000000000000000000000000000000000000000000000! !=====================================! ! Prepare normalization values ! tnorm(2)= CohNorm ! tnorm(3)= NCnorm ! tnorm(4)= CCnorm ! tnorm(5)= obgnorms(ihist,iteration) ! !=====================================! !------------------------------------------------------------------! !()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()! ! Compute Chi^2 array: minchi=99999999999999.0D0 !--> Initialize chisq-minimum variable do inorm = 1, asiz !--> Loop over scalings (DL-01) !=================================================! ! Set temporary floating normalization: ! norm = dble(inorm-1)*istep ! if(ifloat.eq.1) tnorm(3)=norm !--> Float NCDIS ! if(ifloat.eq.2) tnorm(4)=norm !--> Float CCDIS ! if(ifloat.eq.3) tnorm(2)=norm !--> Float CohRho ! !=================================================! chisq(inorm)=0.0D0 do ii=1,nbins(ifloat) !--> Loop over bins ! Subtract Floated MC from data chisqtmp= histvec(1,ii,ifloat) - | histvec(2,ii,ifloat)*tnorm(2) - | histvec(3,ii,ifloat)*tnorm(3) - | histvec(4,ii,ifloat)*tnorm(4) - | histvec(5,ii,ifloat)*tnorm(5) ! Subtract other MC components if(ntyp.ge.6) then do ityp=6,ntyp chisqtmp=chisqtmp-histvec(ityp,ii,ifloat) enddo endif ! Calculate Errors errortmp = dsqrt(histerr(1,ii,ifloat)**2) if(useMCerror.eq.1) then errortmp = dsqrt( errortmp**2 + | (histerr(2,ii,ifloat)*tnorm(2))**2 + | (histerr(3,ii,ifloat)*tnorm(3))**2 + | (histerr(4,ii,ifloat)*tnorm(4))**2 + | (histerr(5,ii,ifloat)*tnorm(5))**2 ) if(ntyp.ge.6) then do ityp=6,ntyp errortmp = dsqrt( errortmp**2 + | histerr(ityp,ii,ifloat)**2 ) enddo endif endif ! Sum Chisq. for each bin if(histvec(1,ii,ifloat).gt.0.0D0) | chisq(inorm) = chisq(inorm) + (chisqtmp**2 / errortmp**2) enddo !--> Loop over bins ! Look for Chi^2 minimum: if(chisq(inorm).lt.minchi) then minchi=chisq(inorm) if(ifloat.eq.1) ncnorms(ihist,iteration,0)=norm if(ifloat.eq.2) ccnorms(ihist,iteration,0)=norm if(ifloat.eq.3)cohnorms(ihist,iteration,0)=norm endif enddo !--> Loop over normalizations !()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()! !------------------------------------------------------------------! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! Find deviation of the Chi^2 Minimum: sigma=dsqrt(minchi/dble(nbins(ifloat)-nzerobins(ifloat))) MaxNorm =-999999999.0D0 MinNorm = 999999999.0D0 do inorm=1,asiz norm = dble(inorm-1)*istep if(dabs((chisq(inorm)-minchi)).lt.sigma) then if(norm.lt.MinNorm) MinNorm=norm if(norm.gt.MaxNorm) MaxNorm=norm endif enddo !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ![][][][][][][][][][][][][][][][][][][][][][][][][]! ! Set Output Values: if(ifloat.eq.1) then ncnorms(ihist,iteration,1) = MaxNorm ncnorms(ihist,iteration,2) = MinNorm ncnorms(ihist,iteration,3) = minchi ncnorms(ihist,iteration,4) = | dble(nbins(ifloat)-nzerobins(ifloat)) elseif(ifloat.eq.2) then ccnorms(ihist,iteration,1) = MaxNorm ccnorms(ihist,iteration,2) = MinNorm ccnorms(ihist,iteration,3) = minchi ccnorms(ihist,iteration,4) = | dble(nbins(ifloat)-nzerobins(ifloat)) elseif(ifloat.eq.3) then cohnorms(ihist,iteration,1) = MaxNorm cohnorms(ihist,iteration,2) = MinNorm cohnorms(ihist,iteration,3) = minchi cohnorms(ihist,iteration,4) = | dble(nbins(ifloat)-nzerobins(ifloat)) endif ![][][][][][][][][][][][][][][][][][][][][][][][][]! !000000000000000000000000000000000000000000000000000000000000000 ! Make output Histogram: ! Fill real array: do ii=1,asiz chisqr(ii)=real(chisq(ii)) enddo ! Set Filename length=LEN_TRIM(outdir(ihist)) if (ifloat.eq.1) | hfile = 'output/'//outdir(ihist)(1:length)//'/chisq_nc.h' if (ifloat.eq.2) | hfile = 'output/'//outdir(ihist)(1:length)//'/chisq_cc.h' if (ifloat.eq.3) | hfile = 'output/'//outdir(ihist)(1:length)//'/chisq_coh.h' ! Open output histogram file: call hropen(60,'OUTPUT',hfile,'N',1024,istat) call hcdir('//OUTPUT',' ') ! Book Histogram: if(ifloat.eq.1) | call hbook1(101,'Chisq. For NCDIS',asiz,0.0,3.0,0.0) if(ifloat.eq.2) | call hbook1(101,'Chisq. For CCDIS',asiz,0.0,3.0,0.0) if(ifloat.eq.3) | call hbook1(101,'Chisq. For CohRho0',asiz,0.0,3.0,0.0) ! Pack Histogram call hpak(101,chisqr) ! Close Output Histograms call HCDIR('//OUTPUT',' ') call hrout(0,icycle, ' ') call hrend('OUTPUT') call hdelet(101) close(60) !000000000000000000000000000000000000000000000000000000000000000 enddo !--> Loop over floated components !|||||||||||||||||||||||||||||||||||||||||||||||||||! !+++++++++++++++++++++++++++++++++++++++++++++++++++! ! ! enddo !--> Loop over ihist ! ! ! !+++++++++++++++++++++++++++++++++++++++++++++++++++! !|||||||||||||||||||||||||||||||||||||||||||||||||||! return end