** Unified Analysis Signal Extractor, as used in PRL for pure D2O. ** Usage: sigex < sigex_cards ** Origin Scott Oser, 25/06/01 ** Modified to run at Carleton by Ferenc and Richard, Feb 2003 ** Remove RSP data energy drift correction when using pseudo MC data ** To use this program for constrained B8 CC shape, set num_cc_bins to 1, ** and adjust number of PDFs. See the accompanying include file * sigex.inc * Log: Modified by AlainB to get Ntuple for the signal MC (May 2003) * Modified by AlainB to get data in ASCI file (Jan 2004) * * PROGRAM UNIFIED_EXTRACTOR implicit none INCLUDE 'sigex.inc' external fcn integer i,j,k character*20 cmdarg,cmdarg1 used_dnbins = 1 used_pdfs = 3 call getarg(1,cmdarg) if ((cmdarg .eq. '') .or. (cmdarg .eq. 'help')) then write (*,*) +'syntax: program [settings] ' stop endif call getarg(2,cmdarg) if (cmdarg .eq. '') then cmdarg='settings_mc' endif write (*,*) 'Reading settings from file ',cmdarg ** Open the file 'settings', which contains a table of response ** function and low energy background settings open(93,file=cmdarg, form='formatted') do i=1,20 read(93, *) j, settings(i,1),k,settings(i,2) if (used_dnbins .eq. 1) then settings(i,1)=settings(i,1)+settings(i,2) settings(i,2)=0 if (i .gt. 5) settings(i,1)=settings(i,1)/2 endif write (*,*) i, settings(i,1), settings(i,2) enddo close(93) ** Maximize the likelihood function with MINUIT ** Uses runcards for inputs to signal extraction call MINUIT (fcn,0) ** Call one last time to calculate chisquared values call fcn (-1) write (*,*) 'End' STOP END ** The fcn routine calculates and returns the likelihood ** for MINUIT. Upon initialization, it also reads in the ** data and PDFs. SUBROUTINE FCN( NPAR,GR,F,X,IFLAG,DPFN ) implicit none INCLUDE 'sigex.inc' real*8 x,f,gr,dpfn,saved_x integer npar,iflag real*8 expected,sum ** Define which energy estimator you want to use. ** Use either RSP or NHIT. Comment out the appropriate ** line, depending on which ntuple definition you're using. ** All you're doing here is telling the program which ntuple ** entry contains the energy variable you want to use for ** the signal extraction. integer ENERGY_EST C parameter (ENERGY_EST = 48) * For ntuple 310, entry 81 is drift-corrected RSP energy parameter (ENERGY_EST = 81) character*20 cmdarg character*128 ntp_fname real*8 totalchisq,denom ** This array stores the binned good events real*8 events dimension events(num_energybins,num_r3bins,num_cstsunbins, + num_dnbins) ** This array stores the PDF shapes themselves real*8 pdfs dimension pdfs(num_pdfs,num_energybins,num_r3bins,num_cstsunbins, + num_dnbins) real*8 bkgd_pdf dimension bkgd_pdf(num_energybins,num_r3bins,num_cstsunbins, + num_dnbins) real*8 pd_background, num_events(num_pdfs) real*8 conv(num_pdfs) ** This arrays are used by MINUIT. x(num_pdfs) contains ** the current best guesses for the number of events for each PDF. ** We save the value of this array between calls, so it will be ** ready at the final call we make to calculate chisquareds. dimension gr(num_pdfs),x(num_pdfs) dimension saved_x(num_pdfs) ** In these arrays we store the radial, energy, and angular profiles ** of the data real*8 r3profile(num_r3bins) real*8 energyprofile(num_energybins) real*8 cstsunprofile(num_cstsunbins) integer firstpass data firstpass /1/ integer i,j,k integer nbin,cbin,rbin,zbin integer cut * Number of actual events of each type, for MC data sets integer actual(num_pdfs) ** Some necessary and obscure HBOOK stuff to read ntuples integer nwpawc,ntuple,ierr, h parameter (nwpawc=2000000) integer record_length /1024/ common/pawc/h(nwpawc) integer ntuple_size ** Use ntuple definition 210 C parameter (ntuple=210) C parameter (ntuple_size = 79) ** Use ntuple definition 310 parameter (ntuple=310) parameter (ntuple_size = 81) integer nevent,ntes,ntnc,ntcc integer good,totalgood real circ real*8 energy,r3,cstsun,npmt,zenith real*8 itr, iretc, thetaij, ang2d, angphi real data_ntp(ntuple_size) real*8 rsp_drift,alt_rsp_drift_1,alt_rsp_drift_2 real azimuth * Ntuple entry definitions for Ntuple 210 and 310 integer TIME_JDY /9/ integer TIME_S /10/ integer NHITS /15/ integer XFP /27/ integer YFP /28/ integer ZFP /29/ integer COSTSUN /36/ integer ITR_WORD /46/ integer THETAIJ_WORD /37/ integer A2DPRB_WORD /34/ integer PHIPRB_WORD /35/ c integer XFG /40/ c integer YFG /41/ c integer ZFG /42/ integer RDAMN /23/ integer DATASET /6/ integer RUNNO /3/ integer DRZSUN /20/ integer GTID /5/ integer INT_PPP /71/ save firstpass, events, actual,pdfs,bkgd_pdf ** Do the following upon the first call to this routine if (firstpass .eq. 1) then ** Zero the event histogram, and the profile histograms do nbin = 1,num_energybins do rbin = 1,num_r3bins do cbin = 1,num_cstsunbins do zbin = 1,used_dnbins events (nbin,rbin,cbin,zbin) = 0 r3profile(rbin) = 0 energyprofile(nbin) = 0 cstsunprofile(cbin) = 0 bkgd_pdf(nbin,rbin,cbin,zbin) = 0 enddo enddo enddo enddo do i=1,used_pdfs actual(i) = 0 enddo ** Reset the counter for the number of good events totalgood = 0 call hlimit(nwpawc) * Get filename from command line CALL getarg(1, ntp_fname) call hropen (1,'exam2', ntp_fname, + 'p',record_length,j) write(*,*) 'Reading file ',ntp_fname * Initialize ntuple reading, then get number of events in ntuple call hgnpar(ntuple,'test') call hnoent(ntuple,nevent) write (*,*) 'Reading data ntuple events =', nevent * Loop over number of events in ntuple * For MC events we can count the number of ES,NC,CC events ntes=0 ntnc=0 ntcc=0 do k=1,nevent call hgnf (ntuple,k,data_ntp,ierr) if(data_ntp(INT_PPP).eq.110)ntes=ntes+1 if(data_ntp(INT_PPP).eq.111)ntnc=ntnc+1 if(data_ntp(INT_PPP).eq.112)ntcc=ntcc+1 enddo write (*,*) 'Total number of CC, ES, NC events =', ntcc,ntes,ntnc ntes=0 ntnc=0 ntcc=0 *** debug alainb write(*,*) ' ' write(*,*) ' energy radius cthsun' do i=1,nevent * This reads event i from ntuple ntuple into array data_ntp call hgnf (ntuple,i,data_ntp,ierr) * Assume this is a good event until proven otherwise good = 1 energy = data_ntp(ENERGY_EST) * If using RSP, apply energy drift to the data! This is the * piecewise linear drift model from Mark Boulay. C if (ENERGY_EST .ne. NHITS) then C energy=energy / rsp_drift(data_ntp(TIME_JDY)) C endif npmt = data_ntp(NHITS) * Calculate r**3 from the path fitter. r3 = (data_ntp(XFP)**2. + data_ntp(YFP)**2. + + data_ntp(ZFP)**2.)**1.5 * Initialize other interesting variables cstsun = data_ntp(COSTSUN) zenith = data_ntp(DRZSUN) itr = data_ntp(ITR_WORD) iretc = 0 thetaij = data_ntp(THETAIJ_WORD) ang2d = data_ntp(A2DPRB_WORD) angphi = data_ntp(PHIPRB_WORD) * Use grid fitter instead, disabling path fitter FOMs c r3 = (data_ntp(XFG)**2. + data_ntp(YFG)**2. c + + data_ntp(ZFG)**2.)**1.5 c angphi = 1 c ang2d = 1 * Apply muon follower short cut--turn off good flag if muon follower if (btest(int(data_ntp(RDAMN)),2)) then good = 0 endif if (btest(int(data_ntp(RDAMN)),11)) then good = 0 endif * Apply day/night cuts, with blindness cut if ((day_night .eq. 1) .and. (data_ntp(DRZSUN) .lt. 0.)) then good = 0 endif if ((day_night .eq. -1) .and. (data_ntp(DRZSUN) .gt. 0.)) then good = 0 endif * Determine which energy, radial, and angular bin the event * falls into--answer is returned in nbin, rbin, cbin call pickbin (energy,r3,cstsun, + zenith,nbin,rbin,cbin,zbin) * See if this event is in our signal box * Does event fall in our fiducial region? * If not, toss it out. if ((nbin .lt. 0) .or. (rbin .lt. 0) .or. (zbin .lt. 0) .or. + (cbin .lt. 0)) then good = 0 endif * Apply data/fit quality cuts if (cut(itr,thetaij,ang2d,angphi,npmt) .ne. 0) good=0 * If good event, bin it! if (good .eq. 1) then * Add event to signal extraction histogram and to * profile histograms. events(nbin,rbin,cbin,zbin) = + events(nbin,rbin,cbin,zbin) + 1 r3profile(rbin) = r3profile(rbin) + 1 energyprofile(nbin) = energyprofile(nbin) + 1 cstsunprofile(cbin) = cstsunprofile(cbin) + 1 *** debug alainb: ASCI file for the data write(*,*) float(nbin)-0.5, float(rbin)-0.5, float(cbin)-0.5 *** debug alainb * Increment good event counter totalgood = totalgood + 1 * Count actual number of MC events of each type accepted * For MC events we can count the number of ES,NC,CC events if(data_ntp(INT_PPP).eq.110)ntes=ntes+1 if(data_ntp(INT_PPP).eq.111)ntnc=ntnc+1 if(data_ntp(INT_PPP).eq.112)ntcc=ntcc+1 c write (*,*) 'good',int(data_ntp(3)), int(data_ntp(5)), c + energy,nbin else c write (*,*) 'bad event',energy,r3,cstsun endif enddo write(*,*) ' ' * This closes hbook file call hrend('exam2') close(1) write (*,*) 'End of reading ntuple' write (*,*) 'Number of good CC, ES, NC events =', ntcc,ntes,ntnc write (*,*) 'Energy profile of data:' do i=1,num_energybins write (*,*) i,energyprofile(i) enddo write (*,*) 'R**3 profile of data:' do i=1,num_r3bins write (*,*) i,r3profile(i) enddo write (*,*) 'Angular profile of data:' do i=1,num_cstsunbins write (*,*) i,cstsunprofile(i) enddo write (*,*) write (*,*) 'Total good events in data =',totalgood * Now that data is read in, now create the PDFs call assemble_pdfs (pdfs) call assemble_bkgd_pdf (bkgd_pdf) * End of things to do on first pass firstpass = 0 write (*,*) 'Ntuple Data Read In!' endif * Calling this routine with npar=-1 is a flag for * calculating the chisquared, using the events and pdfs * arrays, and the saved number of events of each signal type. if (npar .eq. -1) then write (*,*) 'Calculating chisquareds!' call calc_chisqs(events,pdfs,bkgd_pdf,saved_x) write (*,*) ' ' return endif * Calculate Extended Max. Likelihood for data points * See Cowan's statistics book for derivation sum = 0. totalchisq = 0 * Loop over the 3-D event histogram. do nbin = 1, num_energybins do rbin = 1, num_r3bins do cbin = 1, num_cstsunbins do zbin = 1, used_dnbins expected = 0 * Calculate expected number of events in this bin do j=1, used_pdfs num_events(j) = x(j) * Use these when fitting for phi_e, phi_mu directly pd_background=0 if (j .ge. num_cc_bins*used_dnbins+3) then pd_background=settings(1,1+mod(j+1,2)) endif if ((used_dnbins .eq. 1) + .and. (j .eq. num_cc_bins+2)) then pd_background=settings(1,1) *** debug alainb: removing 51.7 internal neutron background * write(*,*) 'Come here',settings(1,1) *** debug alainb endif num_events(j) = num_events(j) + pd_background expected = expected + + num_events(j)*pdfs(j,nbin,rbin,cbin,zbin) c write(*,*) 'pdf',j,nbin,rbin,cbin,zbin c + pdfs(j,nbin,rbin,cbin,zbin),num_events(j) enddo expected=expected+bkgd_pdf(nbin,rbin,cbin,zbin) * If you see this message, your PDF is screwed up, and probably * is not normalized correctly. if (expected .lt. 0) then write (*,*) 'PDF error' write (*,*) x(1),x(2),x(3),x(4) write (*,*) pdfs(1,nbin,rbin,cbin,zbin), + nbin,rbin,cbin,zbin write (*,*) pdfs(2,nbin,rbin,cbin,zbin), + nbin,rbin,cbin,zbin write (*,*) pdfs(3,nbin,rbin,cbin,zbin), + nbin,rbin,cbin,zbin write (*,*) pdfs(4,nbin,rbin,cbin,zbin), + nbin,rbin,cbin,zbin stop endif * Increment sum, which is the sum over events for the log likelihood if (expected .ne. 0) then sum = sum + events(nbin,rbin,cbin,zbin)* + log(expected) * Increment a crude chisquared for the signal extraction. * Chisquared is not the appropriate figure of merit in any case, * so this stuff is better ignored. denom = events(nbin,rbin,cbin,zbin) if (denom .lt. 1.) denom = 1 totalchisq = totalchisq + + (events(nbin,rbin,cbin,zbin)-expected)**2/denom else * If the expected number of events in this bin is zero, but you * actually do have an event there, you have a problem. Basically * you have an event at a place where your fit hypothesis says there * can't be events. You had better fix this. If you're limited by * Monte Carlo statistics, run more Monte Carlo. You can always * cheat and add a single event to one PDF for the appropriate bin. if (events(nbin,rbin,cbin,zbin) .ne. 0) then write (*,*) 'Fatal error!' write (*,*) 'Event present where PDFs are zero' write (*,*) nbin,rbin,cbin,zbin stop endif endif enddo enddo enddo enddo * Since this is extended maximum likelihood, subtract off the sum * of the signal amplitudes do j = 1, used_pdfs sum = sum - num_events(j) enddo * Minuit wants to minimize things, so take the negative of the * log likelihood f = -sum do j = 1,used_pdfs saved_x(j) = num_events(j) enddo * If you give a damned about chisquared for signal extraction, * print it here * write (*,*) 'Total chisquared for signal extraction', totalchisq RETURN END ** This routine calculates bins for data ** It also defines the fiducial volume, energy threshold, and binning ** for the signal extraction subroutine pickbin (energy,r3,cstsun,drzsun,nbin,rbin,cbin,zbin) IMPLICIT NONE INCLUDE 'sigex.inc' integer idum common idum real ran1 real*8 energy,r3,cstsun,drzsun integer nbin,rbin,cbin,zbin * Define the upper and lower energy limits * energy_min is the lower threshold * The quantity (energy_max-energy_min)/num_energybins * defines the size of each energy bin. * The final energy bin is then extended in length * out to absolute_cutoff. real*8 energy_min,energy_max,absolute_cutoff * NHIT-based limits c parameter (energy_min = 46.0) c parameter (energy_max = 172.0) c parameter (absolute_cutoff = 180.) * RSP-based limits for NC PRL parameter (energy_min = 5.511) c parameter (energy_max = 9.011) c parameter (absolute_cutoff = 9.011) parameter (energy_max = 16.011) parameter (absolute_cutoff = 20.011) * Values for previous CC PRL c parameter (energy_min = 7.25) c parameter (energy_max = 16.00) c parameter (absolute_cutoff = 20.00) * Define the minimum and maximum radius for fit real*8 r3_min,r3_max parameter (r3_min = 0.) parameter (r3_max = 166.375e6) real*8 cstsun_min,cstsun_max parameter (cstsun_min = -1.) parameter (cstsun_max = 1.) * Calculate which energy bin the event falls in. * This is not the same as which CC spectrum bin! * Return a negative value if event is out of range nbin = int (num_energybins*(energy-energy_min) + /(energy_max-energy_min))+1 if (energy .le. energy_min) nbin = -1 if (energy .ge. energy_max) nbin = num_energybins if (energy .ge. absolute_cutoff) nbin = -1 rbin = int (num_r3bins*1.0*(r3-r3_min)/(r3_max-r3_min))+1 if (r3 .lt. r3_min) rbin = -1 if (r3 .ge. r3_max) rbin = -1 * Calculate which angular bin the event falls in * This commented-out code uses equal angular bins from -1 to +1 cbin = int (num_cstsunbins* + (cstsun-cstsun_min)/(cstsun_max-cstsun_min))+1 * The following puts half the bins from 0.5-1.0 in cos(theta) * This gives better coverage over the ES peak. c if (cstsun .le. 0.5) cbin = int (num_cstsunbins*0.5* c + (cstsun-cstsun_min)/(0.5-cstsun_min))+1 c if (cstsun .gt. 0.5) cbin = num_cstsunbins/2 + c + int (num_cstsunbins*0.5* c + (cstsun-0.5)/(cstsun_max-0.5))+1 if (cstsun .lt. cstsun_min) cbin = -1 if (cstsun .gt. cstsun_max) cbin = -1 if (drzsun .lt. 0) then zbin = 1 else zbin = 2 endif if (zbin .gt. used_dnbins) zbin=used_dnbins return end ** This routine reads in and assembles the PDFs subroutine assemble_pdfs (pdfs) IMPLICIT NONE INCLUDE 'sigex.inc' real*8 pdfs dimension pdfs(num_pdfs,num_energybins,num_r3bins,num_cstsunbins, + num_dnbins) real*8 energy,r3,cstsun, normalize,npmt, runnum, zenith integer cut integer i,j,k,l,m,current_pdf integer record_length /1024/ real*8 energyprof(num_energybins),cc_energyprof(num_energybins) real*8 r3prof(num_r3bins) real*8 cstsunprof(num_cstsunbins) c integer energyprof(num_energybins) c integer r3prof(num_r3bins) c integer cstsunprof(num_cstsunbins) integer nbin,rbin,cbin,zbin character*127 pdf_filename real*8 pee, weight integer idum common idum real ran1 * Number of CC spectrum bins. Set to 1 to use B8 shape constraint * Else set to number of bins you want in the CC spectrum integer cc_bin real sigma, sigmaerr, gasdev real*8 rescale_nhit, rescale_rsp, mc_scale * Required HBOOK variables for reading ntuples integer nwpawc,ntuple,ierr, h parameter (nwpawc=2000000) integer nevent integer ntuple_size integer good real*8 totalgood real circ integer whichfile * Ntuple parameters for Unified Ntuple c parameter (ntuple=210) c parameter (ntuple_size = 79) parameter (ntuple=310) parameter (ntuple_size = 81) * Specify which energy estimator is being used. integer ENERGY_EST parameter (ENERGY_EST = 48) * Ntuple entry definitions for Ntuple 210 integer TIME_JDY /9/ integer TIME_S /10/ integer NHITS /15/ integer XFP /27/ integer YFP /28/ integer ZFP /29/ integer COSTSUN /36/ integer ITR_WORD /46/ integer THETAIJ_WORD /37/ integer A2DPRB_WORD /34/ integer PHIPRB_WORD /35/ c integer XFG /40/ c integer YFG /41/ c integer ZFG /42/ integer RDAMN /23/ integer DATASET /6/ integer RUNNO /3/ integer DRZSUN /20/ integer GTID /5/ integer NUE /68/ integer loop, how_many_times real azimuth *** alainb debug integer tagdim parameter (tagdim=5) REAL TAGVAL(tagdim) CHARACTER*6 TAGSIG(tagdim) DATA TAGSIG /'energy','radius','cthsun','thetij','tagtyp'/ *** alainb debug real*8 itr, iretc, thetaij, ang2d, angphi real data_ntp(ntuple_size) common/pawc/h(nwpawc) * For technical reasons, we require an even number of angular bins.. * This could be changed by rewriting pickbin accordingly. if (int(num_cstsunbins*0.5) .ne. num_cstsunbins*0.5) then write (*,*) 'Please use an even number of cos theta bins' c stop endif if (used_pdfs .lt. num_cc_bins) then write (*,*) 'Error: You are using fewer PDFs than CC bins!' stop endif if (num_cc_bins .eq. 1) then write (*,*) 'Using B8 CC spectrum as constraint' else write (*,*) 'Not constraining CC spectral shape' endif * First, zero all the pdfs do current_pdf = 1,used_pdfs do i = 1, num_energybins do j = 1, num_r3bins do k = 1, num_cstsunbins do l = 1, used_dnbins pdfs(current_pdf,i,j,k,l) = 0 enddo enddo enddo enddo enddo call hlimit(nwpawc) *** debug alainb call hcdir('//PAWC',' ') call hmdir('user','s') call hbook1(1001,'cc ebin', + num_energybins,0.,float(num_energybins),0.) call hbook1(1002,'cc rbin', + num_r3bins,0.,float(num_r3bins),0.) call hbook1(1003,'cc cbin', + num_cstsunbins,0.,float(num_cstsunbins),0.) CALL HBOOKN(1111,'cc ntuple ',TAGDIM,'user',1024, TAGSIG) call hidopt(0,'stat') call hbook1(2001,'es ebin', + num_energybins,0.,float(num_energybins),0.) call hbook1(2002,'es rbin', + num_r3bins,0.,float(num_r3bins),0.) call hbook1(2003,'es cbin', + num_cstsunbins,0.,float(num_cstsunbins),0.) CALL HBOOKN(2222,'es ntuple ',TAGDIM,'user',1024, TAGSIG) call hidopt(0,'stat') call hbook1(3001,'nc ebin', + num_energybins,0.,float(num_energybins),0.) call hbook1(3002,'nc rbin', + num_r3bins,0.,float(num_r3bins),0.) call hbook1(3003,'nc cbin', + num_cstsunbins,0.,float(num_cstsunbins),0.) CALL HBOOKN(3333,'nc ntuple ',TAGDIM,'user',1024, TAGSIG) call hidopt(0,'stat') *** debug alainb * Open the file 'pdflist210', which contains a list of all the * PDF filenames. if (used_dnbins .eq. 1) then open(93,file='pdflist210', form='formatted') else open(93,file='pdflist210_dn', form='formatted') endif * Now loop over pdfs current_pdf = 0 do while (current_pdf .lt. used_pdfs) do m=1,used_dnbins current_pdf = current_pdf + 1 * Zero the profile histograms do i = 1, num_energybins energyprof(i) = 0 enddo do i = 1, num_r3bins r3prof(i) = 0 enddo do i = 1, num_cstsunbins cstsunprof(i) = 0 enddo totalgood = 0 do whichfile=1,3 * The PDFs are assumed to reside in either the directory * Current_NHIT_PDFs or Current_RSP_PDFs c if ((current_pdf .eq. 1) .or. (current_pdf .gt. num_cc_bins)) then read(93, *) pdf_filename if (ENERGY_EST .eq. NHITS) then pdf_filename = 'Oct_NHIT_PDFs/'//pdf_filename else c pdf_filename = '/tape/sno_disk/signal_mc/June2/PDFs/' c + //pdf_filename pdf_filename = 'Oct_RSP_PDFs/'//pdf_filename endif c endif * Open the ntuple file for this pdf call hropen (1,'exam2', pdf_filename, + 'p',record_length,j) call hgnpar(ntuple,'test') call hnoent(ntuple,nevent) write (*,*) 'Reading pdf',pdf_filename do i=1,nevent * This reads event i from ntuple ntuple into array data_ntp call hgnf (ntuple,i,data_ntp,ierr) good = 1 zenith = data_ntp(DRZSUN) * Call initially just to get zbin right call pickbin (0.d0,0.d0,0.d0,zenith,nbin,rbin,cbin,zbin) energy = data_ntp(ENERGY_EST) * If using NHIT, make it a continuous variable first. * This is because NHIT is naturally a discrete variable, but * we may need to rescale it. We do this by making it continuous * by adding a random number from -0.5 to +0.5, rescaling, then * rounding off. if (ENERGY_EST .eq. NHITS) energy = energy - 0.5 + ran1(idum) npmt = data_ntp(NHITS) runnum = data_ntp(RUNNO) c The scale factor used for MC generation if (ENERGY_EST .eq. NHITS) then mc_scale = 0.5834 else mc_scale = 0.5672 endif * Andre's recommended energy drift for NHIT--May 21, 2001 rescale_nhit = (1.01553 - 5.523e-5*(data_ntp(TIME_JDY)-9072)) + * mc_scale * Andre's recommended scaling value for RSP with 4.0186 Unified MC * If rescale_rsp = mc_scale, then no rescaling necessary! rescale_rsp = 0.5672 * Rescale the energy. This is mandatory for the NHIT analysis. * For the RSP analysis, we correct the data, and not the MC. * So for RSP, rescale_rsp should normally equal mc_scale if (ENERGY_EST .eq. NHITS) then energy = energy*rescale_nhit/mc_scale else energy = energy*rescale_rsp/mc_scale endif npmt = npmt*rescale_nhit/mc_scale * Scale energy scale of PDF energy = energy*settings(6,zbin) * If energy resolution smearing enabled, add Gaussian error to energy c sigma = .1432 + .14488*data_ntp(67) c + -.002357*data_ntp(67)*data_ntp(67) sigma = -0.06837+0.3308*sqrt(data_ntp(67)-.511) + +.04253*(data_ntp(67)-.511) c sigmaerr=sigma*sqrt((1.045+.00401*(data_ntp(67)-5.49))**2. - 1.) sigmaerr=sigma*sqrt((1.045+.00401*(energy-5.486))**2. - 1.) if (current_pdf .eq. 3) sigmaerr=0.314 if (settings(7,zbin) .eq. 1) then energy=energy+gasdev(idum)*sigmaerr endif c Discretize energy if using NHIT if (ENERGY_EST .eq. NHITS) energy = int(energy+0.5) npmt = int(npmt+0.5) * smear vertex resolution data_ntp(XFP)=data_ntp(XFP)+settings(10,zbin)*gasdev(idum) data_ntp(YFP)=data_ntp(YFP)+settings(10,zbin)*gasdev(idum) data_ntp(ZFP)=data_ntp(ZFP)+settings(10,zbin)*gasdev(idum) * Use path fitter r3 = (data_ntp(XFP)**2. + data_ntp(YFP)**2. + + data_ntp(ZFP)**2.)**1.5 * Scale in radius r3 = r3/((settings(9,zbin))**3.) cstsun = data_ntp(COSTSUN) itr = data_ntp(ITR_WORD) iretc = 0 thetaij = data_ntp(THETAIJ_WORD) * corrchisq = data_ntp(no_corrchisq) ang2d = data_ntp(A2DPRB_WORD) angphi = data_ntp(PHIPRB_WORD) * let's smear cos(sun) if (ran1(idum) .lt. settings(11,zbin)) then cstsun=2.*(ran1(idum)-0.5) endif * Rotate sun direction by settings(12) degrees cstsun = cstsun*cos(settings(12,zbin)*3.14159/180) + +sin(settings(12,zbin)*3.14159/180)* + sqrt(1-cstsun*cstsun)*sin(6.28*ran1(idum)) * Add in Aksel crosstalk smearing function correction for RSP if (settings(8,zbin) .eq. 1) then energy = energy + .0023*19.1*(energy-5.49)/(19.1-5.49) endif call pickbin (energy,r3,cstsun,zenith,nbin,rbin,cbin,zbin) * See if this event is in our signal box * Does event fall in our fiducial region? if ((nbin .lt. 0) .or. (rbin .lt. 0) .or. (zbin .lt. 0) .or. + (cbin .lt. 0)) then good = 0 endif * Split Day/night PDFs if (zbin .ne. m) good=0 * ITR.le. 0 indicates bad fit--throw out such events if (itr .le. 0) good=0 * Apply data/fit quality cuts. For the Unified Analysis, HLC * were not applied to the Monte Carlo. c if (cut(itr,thetaij,ang2d,angphi,npmt) .ne. 0) good=0 * Eliminate high radon runs c if ((runnum .ge. 10000) .and. (runnum .le. 10749)) good=0 c if ((runnum .ge. 15065) .and. (runnum .le. 15270)) good=0 * Use only runs through 12168 c if (runnum .gt. 12168) good=0 * Use only runs through 14685 c if (runnum .gt. 14685) good=0 * Cut on east/west c if (azimuth(data_ntp(18),data_ntp(19)) .lt. 3.14159265) c + good=0 * Apply the RO cut c if (circ(data_ntp(TIME_JDY),data_ntp(TIME_S),2.) .ne. 0.) c + good=0 * Apply day/night cuts, with blindness cut if ((day_night .eq. 1) .and. (data_ntp(DRZSUN) .lt. 0.)) then good = 0 endif if ((day_night .eq. -1) .and. (data_ntp(DRZSUN) .gt. 0.)) then good = 0 endif C If good event, bin it! if (good .eq. 1) then weight = 1 *** debug alainb TAGVAL(1) = energy TAGVAL(2) = r3 TAGVAL(3) = cstsun TAGVAL(4) = thetaij if (current_pdf .eq. 1) then call hfill(1001,float(nbin)-0.5,0.,1.) call hfill(1002,float(rbin)-0.5,0.,1.) call hfill(1003,float(cbin)-0.5,0.,1.) TAGVAL(5) = 1. CALL HFN(1111,TAGVAL) endif if (current_pdf .eq. 2) then call hfill(2001,float(nbin)-0.5,0.,1.) call hfill(2002,float(rbin)-0.5,0.,1.) call hfill(2003,float(cbin)-0.5,0.,1.) TAGVAL(5) = 2. CALL HFN(2222,TAGVAL) endif if (current_pdf .eq. 3) then call hfill(3001,float(nbin)-0.5,0.,1.) call hfill(3002,float(rbin)-0.5,0.,1.) call hfill(3003,float(cbin)-0.5,0.,1.) TAGVAL(5) = 3. CALL HFN(3333,TAGVAL) endif *** debug alainb energyprof(nbin) = energyprof(nbin) +weight r3prof(rbin) = r3prof(rbin) +weight cstsunprof(cbin) = cstsunprof(cbin) +weight totalgood = totalgood +weight endif enddo ! endloop on nevent c This closes hbook file call hrend('exam2') close(1) enddo ! end loop on whichfile write (*,*) 'Total good events, pdf', current_pdf, '=',totalgood pdf_count(current_pdf) = totalgood if (current_pdf .eq. 1) then do i=1,num_energybins cc_energyprof(i) = energyprof(i) enddo endif how_many_times = 1 if (current_pdf .eq. 1) how_many_times = num_cc_bins do loop=1,how_many_times * Mask off the energy dependence of the CC pdfs if doing spectrum * It's important to do this before renormalizing PDF. * The routine cc_bin(i) defines how the different energy bins are * grouped into CC spectral bins. * The individual CC spectrum PDFs are generated by taking the regular * constrained CC PDF, and replacing the energy component with a * new energy shape--flat inside the spectral bin, and zero outside! if ((num_cc_bins .gt. 1).and.(current_pdf .le. num_cc_bins)) then pdf_count(current_pdf) = 0. do i = 1, num_energybins if (cc_bin(i) .eq. current_pdf) then pdf_count(current_pdf) = + pdf_count(current_pdf)+cc_energyprof(i) energyprof(i) = 1 else energyprof(i) = 0 endif enddo endif * The PDF is assumed to factorize into energy, radial, and angular * components. So calculate the 3-D PDF shape by multiplying, then * normalize to unit probability. normalize = 0 do i = 1, num_energybins do j = 1, num_r3bins do k = 1, num_cstsunbins do l = 1, used_dnbins pdfs(current_pdf,i,j,k,l) = (energyprof(i)) + * (r3prof(j))*(cstsunprof(k)) + /(totalgood**3.) if (l .ne. m) pdfs(current_pdf,i,j,k,l)=0 normalize = normalize + pdfs(current_pdf,i,j,k,l) enddo enddo enddo enddo write (*,*) 'Normalize =', normalize * Write out the PDF shape. Note that this assumes that * num_energybins .ge. num_r3bins or num_cstsunbins. * Note that if num_energybins .gt. the other two, then the * last few entries of those columns will be junk. do i = 1, num_energybins write (*,*) current_pdf,i,energyprof(i),r3prof(i),cstsunprof(i) enddo * Finish normalizing PDFs. do i = 1, num_energybins do j = 1, num_r3bins do k = 1, num_cstsunbins do l = 1, used_dnbins pdfs(current_pdf,i,j,k,l) = + pdfs(current_pdf,i,j,k,l) / normalize enddo enddo enddo enddo if (current_pdf .lt. num_cc_bins) current_pdf = current_pdf+1 * End "loop" loop on how many times to produce the PDF! enddo write (*,*) 'At end of this pass, current_pdf = ', current_pdf C End loop over dnbins enddo C End loop over pdfs (while loop) enddo close(93) *** debug alainb call hrput(0,'ryh.hist','n') *** debug alainb return end * Apply High Level Cuts. Return 1 if event is to be cut, * else return 0 for good events. function cut (itr, thetaij, ang2d, angphi, npmt) IMPLICIT NONE integer cut real*8 itr, thetaij, ang2d, angphi,npmt cut = 0 if (itr .le. 0.55) cut = 1 if (thetaij .ge. 1.45) cut = 1 if (thetaij .le. 0.75) cut = 1 if (itr .le. 0.0) cut = 1 if (ang2d .le. 4.e6/npmt**6) cut = 1 if (angphi .le. 1.e-4) cut = 1 return end * This routine determines which CC spectrum bin an event belongs to, * and puts it there. It returns -1 if event is out of range * Events are classified based on what energy bin they fall into. * For the NC Analysis PRL, the first 26 energy bins were grouped by * pairs into 13 CC spectrum bins, and the final 14 energy bins gave the * 14th spectrum bin. There's nothing sacred about this, however. function cc_bin (nbin) IMPLICIT NONE INCLUDE 'sigex.inc' integer cc_bin integer nbin * Make 11 bins total--the last one is wider than the others c if (nbin .le. 26) cc_bin = int(0.51+nbin*0.5) c if ((nbin .ge. 27) .and. (nbin .le. 40)) cc_bin = 14 if (nbin .le. 2*(num_cc_bins-1)) cc_bin = int(0.51+nbin*0.5) if ((nbin .gt. 2*(num_cc_bins-1)) .and. + (nbin .le. num_energybins)) cc_bin = num_cc_bins * These lines were used for first PRL c if (nbin .le. 20) cc_bin = int(0.51+nbin*0.5) c if ((nbin .ge. 21) .and. (nbin .le. 34)) cc_bin = 11 if ((nbin .lt. 1.) .or. (nbin .gt. num_energybins)) cc_bin=-1 return end * This routine calculates a very crude chisquared between the data * and the Monte Carlo for the radial, energy, and angular distributions. * It's really crude--no account is taken of statistical correlations * between the numbers of extracted events, or uncertainties on the PDFs * themselves. Scott has written a better standalone program to do * these calculations correctly. subroutine calc_chisqs (events, pdfs, bkgd_pdf, x) IMPLICIT NONE INCLUDE 'sigex.inc' real*8 events dimension events(num_energybins,num_r3bins,num_cstsunbins) real*8 pdfs dimension pdfs(num_pdfs,num_energybins,num_r3bins,num_cstsunbins) real*8 bkgd_pdf dimension bkgd_pdf(num_energybins,num_r3bins,num_cstsunbins) real*8 errors(num_pdfs),eplus,eminus real*8 covariance(num_pdfs,num_pdfs) real*8 x(num_pdfs) real*8 chisq,expected,chisqq integer i,j,k,l real*8 energyprofile_data(num_energybins) real*8 energyprofile_mc(num_energybins,num_pdfs+1) real*8 r3profile_data(num_r3bins) real*8 r3profile_mc(num_r3bins,num_pdfs+1) real*8 cstsunprofile_data(num_cstsunbins) real*8 cstsunprofile_mc(num_cstsunbins,num_pdfs+1) real*8 sigma_data,sigma_pred,sigma real*8 bnd1,bnd2,ivar real*8 total_cc_error, total_cc_events,total_cc_prediction real*8 covar_cc_nc,covar_cc_es character*8 lab real*8 conv(num_pdfs) real*8 a,da,delta,night,day,dnight,dday,cov * Conversion constants - calculated in March 02 do i=1,used_pdfs if (i .le. used_dnbins*(num_cc_bins)) conv(i) = 50.1835 if (i .ge. used_dnbins*num_cc_bins+1) conv(i) = 52.4976 if (i .ge. used_dnbins*num_cc_bins+1+used_dnbins) + conv(i) = 51.3865 * This corrects for different eccentricity, livetime - March 02 if (used_dnbins .gt. 1) then if (mod(i,2) .eq. 1) then conv(i) = conv(i) * 0.99469 else conv(i) = conv(i) * 1.00748 endif endif enddo write (*,*) 'Flux conversion factors for SNO pure D2O data' do i=1,used_pdfs write (*,*) i, conv(i) enddo do i=1,used_pdfs call mnpout (i, lab, x(i), sigma,bnd1,bnd2,ivar) enddo call mnemat (covariance, num_pdfs) do i=1,used_pdfs errors(i) = sqrt(covariance(i,i)) enddo c do l=1,used_pdfs c write (*,*) 'Signal', l,x(l), errors(l) c enddo write (*,5) x(1),errors(1),x(2),errors(2), + x(3), errors(3), x(4), errors(4) 5 format (' Final Fit ',8f7.2) total_cc_error = 0 total_cc_events = 0 total_cc_prediction = 0 do i=1,num_cc_bins total_cc_events = total_cc_events+x(i) total_cc_prediction = total_cc_prediction+pdf_count(i) do j=1,num_cc_bins total_cc_error = total_cc_error+covariance(i,j) enddo enddo total_cc_error = sqrt(total_cc_error) write (*,*) 'Total CC error (sum of all bins): ', total_cc_events, + total_cc_error covar_cc_nc = 0 do i=1,num_cc_bins covar_cc_nc = covar_cc_nc + covariance(i,num_cc_bins+2) enddo covar_cc_nc = covar_cc_nc/total_cc_error + /errors(num_cc_bins+2) covar_cc_es = 0 do i=1,num_cc_bins covar_cc_es = covar_cc_es + covariance(i,num_cc_bins+1) enddo covar_cc_es = covar_cc_es/total_cc_error + /errors(num_cc_bins+1) write (*,*) 'CC-NC correlation: ', covar_cc_nc write (*,*) 'CC-ES correlation: ', covar_cc_es do l=1,used_pdfs+1 do i=1,num_energybins energyprofile_data(i) = 0 energyprofile_mc(i,l) = 0 enddo do i=1,num_r3bins r3profile_data(i) = 0 r3profile_mc(i,l) = 0 enddo do i=1,num_cstsunbins cstsunprofile_data(i) = 0 cstsunprofile_mc(i,l) = 0 enddo enddo do i = 1,num_energybins do j = 1,num_r3bins do k = 1,num_cstsunbins energyprofile_data(i) = energyprofile_data(i) + +events(i,j,k) r3profile_data(j) = r3profile_data(j) + +events(i,j,k) cstsunprofile_data(k) = cstsunprofile_data(k) + +events(i,j,k) enddo enddo enddo do l=1,used_pdfs do i = 1,num_energybins do j = 1,num_r3bins do k = 1,num_cstsunbins energyprofile_mc(i,l) = energyprofile_mc(i,l) + +pdfs(l,i,j,k) r3profile_mc(j,l) = r3profile_mc(j,l) + +pdfs(l,i,j,k) cstsunprofile_mc(k,l) = cstsunprofile_mc(k,l) + +pdfs(l,i,j,k) enddo enddo enddo enddo l=used_pdfs+1 do i = 1,num_energybins do j = 1,num_r3bins do k = 1,num_cstsunbins energyprofile_mc(i,l) = energyprofile_mc(i,l) + +bkgd_pdf(i,j,k) r3profile_mc(j,l) = r3profile_mc(j,l) + +bkgd_pdf(i,j,k) cstsunprofile_mc(k,l) = cstsunprofile_mc(k,l) + +bkgd_pdf(i,j,k) enddo enddo enddo * KLUDGE x(3) = x(3) + settings(1,1) chisq = 0 write (*,*) +' Variable Data Fit CC ES NC bkgd nsig' do i=1,num_energybins expected = 0 do l=1,used_pdfs expected = expected + energyprofile_mc(i,l)*x(l) enddo * Add background contribution as well expected = expected + energyprofile_mc(i,used_pdfs+1) sigma_data = sqrt(energyprofile_data(i)) if (sigma_data .eq. 0) sigma_data = 1 sigma_pred = 0 do j=1,used_pdfs do k=1,used_pdfs sigma_pred = sigma_pred + + covariance(j,k)*energyprofile_mc(i,j)*energyprofile_mc(i,k) enddo enddo sigma = sigma_data c Only total chisq for bins with more than 5 events chisqq = (energyprofile_data(i)-expected)/sigma if (energyprofile_data(i) .gt. 5.)chisq = chisq + chisqq**2 write (*,11) 'Energy',energyprofile_data(i),expected, + energyprofile_mc(i,1)*x(1),energyprofile_mc(i,2)*x(2), + energyprofile_mc(i,3)*x(3),energyprofile_mc(i,4),chisqq 11 format(A10,7F8.1) enddo chisq = chisq write (*,*) 'Energy chisq(>5 evs)/nbins=', + chisq,num_energybins chisq = 0 write (*,*) +' Variable Data Fit CC ES NC bkgd nsig' do i=1,num_r3bins expected = 0 do l=1,used_pdfs expected = expected + r3profile_mc(i,l)*x(l) enddo * Add background contribution as well expected = expected + r3profile_mc(i,used_pdfs+1) sigma_data = sqrt(r3profile_data(i)) if (sigma_data .eq. 0) sigma_data = 1 sigma_pred = 0 do j=1,used_pdfs do k=1,used_pdfs sigma_pred = sigma_pred + + covariance(j,k)*r3profile_mc(i,j)*r3profile_mc(i,k) enddo enddo sigma = sigma_data chisqq = (r3profile_data(i)-expected)/sigma chisq = chisq + chisqq**2 write (*,11) 'R3',r3profile_data(i),expected, + r3profile_mc(i,1)*x(1),r3profile_mc(i,2)*x(2), + r3profile_mc(i,3)*x(3),r3profile_mc(i,4),chisqq enddo chisq = chisq write (*,*) 'R3 chisq/nbins=', chisq,num_r3bins chisq = 0 write (*,*) +' Variable Data Fit CC ES NC bkgd nsig' do i=1,num_cstsunbins expected = 0 do l=1,used_pdfs expected = expected + cstsunprofile_mc(i,l)*x(l) enddo * Add background contribution as well expected = expected + cstsunprofile_mc(i,used_pdfs+1) sigma_data = sqrt(cstsunprofile_data(i)) if (sigma_data .eq. 0) sigma_data = 1 sigma_pred = 0 do j=1,used_pdfs do k=1,used_pdfs sigma_pred = sigma_pred + + covariance(j,k)*cstsunprofile_mc(i,j)*cstsunprofile_mc(i,k) enddo enddo sigma = sigma_data chisqq = (cstsunprofile_data(i)-expected)/sigma chisq = chisq + chisqq**2 write (*,11) 'Cstsun',cstsunprofile_data(i),expected, + cstsunprofile_mc(i,1)*x(1),cstsunprofile_mc(i,2)*x(2), + cstsunprofile_mc(i,3)*x(3),cstsunprofile_mc(i,4),chisqq enddo chisq = chisq write (*,*) 'Cstsun chisq/nbins=', chisq,num_cstsunbins * KLUDGE x(3) = x(3) - settings(1,1) do l=1,used_pdfs write (*,6) l,x(l), errors(l), pdf_count(l), + x(l)/pdf_count(l)*conv(l), + errors(l)/pdf_count(l)*conv(l) enddo 6 format ('Signal ',i2,' ',5f12.4) if (used_dnbins .eq. 2) then do l=1,used_pdfs/2 night = x(l*2-1)/pdf_count(l*2-1)*conv(l*2-1) day = x(l*2)/pdf_count(l*2)*conv(l*2) dnight = errors(l*2-1)/pdf_count(l*2-1)*conv(l*2-1) dday = errors(l*2)/pdf_count(l*2)*conv(l*2) a = 2*(night-day)/(night+day) cov = covariance(l*2-1,l*2)/errors(l*2-1)/errors(l*2) delta = dnight**2 + dday**2 - 2*cov*dnight*dday delta = sqrt(delta) da = a*delta/(night-day) write (*,*) 'A:', a, da, night-day, delta, a/da enddo endif return end * The RSP energy drift model function rsp_drift (jdy) real*8 rsp_drift real jdy * First PRL's drift model c if (jdy .lt. 9356) then c rsp_drift = 1.68708 -7.34376e-5*jdy c else c rsp_drift = 1. c endif * New Feb 2002 drift model if (jdy .lt. 9356) then rsp_drift = 1.58719 - 6.28308e-5*jdy else rsp_drift = 0.999345 - 9.12406e-6*(jdy-9356) endif rsp_drift = rsp_drift*1.00506 return end * Alternative RSP energy drift function alt_rsp_drift_1 (jdy) real*8 alt_rsp_drift_1 real jdy c if (jdy .le. 9213) alt_rsp_drift_1 = 1.02233 c if ((jdy .gt. 9213) .and. (jdy .le. 9330)) then c alt_rsp_drift_1 = 1.0066-(jdy-9213.)*1.14515e-4 c endif c if ((jdy .gt. 9330) .and. (jdy .le. 9397)) then c alt_rsp_drift_1 = .99319 c endif c if (jdy .gt. 9397) alt_rsp_drift_1 = 1.0049 if (jdy .le. 9100) alt_rsp_drift_1 = 1.0184 if ((jdy .gt. 9100) .and. (jdy .le. 9212)) then alt_rsp_drift_1 = 1.0184-(jdy-9100.)*8.9286e-5 endif if ((jdy .gt. 9212) .and. (jdy .le. 9394)) then alt_rsp_drift_1 = 1.0084-(jdy-9212.)*9.45055e-5 endif if (jdy .ge. 9394) alt_rsp_drift_1 = 1.0041 if (jdy .ge. 9577) alt_rsp_drift_1 = 0.9939 alt_rsp_drift_1 = alt_rsp_drift_1*1.00506 return end * Alternative RSP energy drift function alt_rsp_drift_2 (jdy) real*8 alt_rsp_drift_2 real jdy c if (jdy .le. 9072) alt_rsp_drift_2 = 1.0156 c if ((jdy .gt. 9072) .and. (jdy .le. 9213)) then c alt_rsp_drift_2 = 1.0156-(jdy-9072.)*6.3830e-5 c endif c if ((jdy .gt. 9213) .and. (jdy .le. 9397)) then c alt_rsp_drift_2 = 1.0126-(jdy-9213.)*4.18478e-5 c endif c if (jdy .gt. 9397) alt_rsp_drift_2 = 0.99319 if (jdy .le. 9100) alt_rsp_drift_2 = 1.0122 if ((jdy .gt. 9100) .and. (jdy .le. 9212)) then alt_rsp_drift_2 = 1.0122-(jdy-9100.)*3.39286e-5 endif if ((jdy .gt. 9212) .and. (jdy .le. 9394)) then alt_rsp_drift_2 = 1.0084-(jdy-9212.)*2.36264e-5 endif if (jdy .ge. 9394) alt_rsp_drift_2 = 0.9950 if (jdy .ge. 9577) alt_rsp_drift_2 = 0.9999 alt_rsp_drift_2 = alt_rsp_drift_2*1.00506 return end function pee (nue) real*8 pee real nue * This is approximately tan^2 theta = 0.3, dm^2 = 2e-5 c pee = .256 + .074*(nue-6)/14. c pee = nue pee = 0.45 - 0.15*((nue/10.)-1) return end subroutine assemble_bkgd_pdf (pdf) IMPLICIT NONE INCLUDE 'sigex.inc' real*8 pdf, component dimension pdf(num_energybins,num_r3bins,num_cstsunbins, + num_dnbins) dimension component(num_energybins,num_r3bins,num_cstsunbins, + num_dnbins) real*8 energyprofile(num_energybins) real*8 r3profile(num_r3bins) real*8 dnprofile(num_dnbins) real*8 e,r3 real*8 normalization real*8 mu,sigma real*8 external_neutrons, d2o_bgs, pmt_bgs, av_bgs real*8 fidvol /0.77d0/ real*8 ethresh /5.5d0/ integer nbin, rbin,cbin, zbin * Required HBOOK variables for reading ntuples integer nwpawc parameter (nwpawc=2000000) *** debug alainb call hlimit(nwpawc) call hcdir('//PAWC',' ') call hmdir('user','s') call hbook1(4001,'bgs ebin', + num_energybins,0.,float(num_energybins),0.) call hbook1(4002,'bgs rbin', + num_r3bins,0.,float(num_r3bins),0.) call hbook1(4003,'bgs cbin', + num_cstsunbins,0.,float(num_cstsunbins),0.) call hidopt(0,'stat') *** debug alainb dnprofile(1) = settings(2,1) dnprofile(2) = settings(2,2) external_neutrons = settings(2,1)+settings(2,2) d2o_bgs = settings(3,1)+settings(3,2) pmt_bgs = settings(4,1)+settings(4,2) av_bgs = settings(5,1)+settings(5,2) write (*,*) 'Assembling background PDF' write (*,*) 'Subtracting',settings(1,1)+settings(1,2), + ' internal p-d neutrons' write (*,*) 'Subtracting',external_neutrons,' external neutrons' write (*,*) 'Subtracting',d2o_bgs,' D2O beta-gammas' write (*,*) 'Subtracting',pmt_bgs,' PMT beta-gammas' write (*,*) 'Subtracting',av_bgs,' AV beta-gammas' write (*,*) ' ' * Zero the background PDF array do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins pdf(nbin,rbin,cbin,zbin)= 0. enddo enddo enddo enddo * Calculate the external neutron component dnprofile(1) = settings(2,1) dnprofile(2) = settings(2,2) mu = 5.5 sigma = 1.03 do nbin=1,num_energybins * This is hard-wired to our analysis threshold e = ethresh-.125+nbin*0.25-mu energyprofile(nbin)=exp(-e*e/2/sigma/sigma) enddo do rbin=1,num_r3bins r3 = (1.0*rbin)/num_r3bins*fidvol r3profile(rbin) = r3 enddo normalization=0. do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins component(nbin,rbin,cbin,zbin)= + energyprofile(nbin)*r3profile(rbin)*dnprofile(zbin) normalization = normalization + +component(nbin,rbin,cbin,zbin) enddo enddo enddo enddo do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins if (normalization .gt. 0) + pdf(nbin,rbin,cbin,zbin) = pdf(nbin,rbin,cbin,zbin) + +component(nbin,rbin,cbin,zbin) + *external_neutrons/normalization enddo enddo enddo enddo * Now calculate the D2O beta-gamma component * Energy spectrum of beta-gammas is tail of Gaussian, according to Josh dnprofile(1) = settings(3,1) dnprofile(2) = settings(3,2) mu = 2.588 sigma = 0.7828 do nbin=1,num_energybins e = ethresh-.125+nbin*0.25-mu energyprofile(nbin)=exp(-e*e/2/sigma/sigma) enddo * Flat R3 distribution for D2O beta-gammas do rbin=1,num_r3bins r3 = (1.0*rbin)/num_r3bins*fidvol r3profile(rbin) = 1 enddo normalization=0. do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins component(nbin,rbin,cbin,zbin)= + energyprofile(nbin)*r3profile(rbin)*dnprofile(zbin) normalization = normalization + +component(nbin,rbin,cbin,zbin) enddo enddo enddo enddo do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins if (normalization .gt. 0) + pdf(nbin,rbin,cbin,zbin) = pdf(nbin,rbin,cbin,zbin) + +component(nbin,rbin,cbin,zbin) + *d2o_bgs/normalization enddo enddo enddo enddo * Now calculate the PMT beta-gamma component * Energy spectrum of PMT beta-gammas from Vadim's Gaussian fit dnprofile(1) = settings(4,1) dnprofile(2) = settings(4,2) mu = 1.416 sigma = 0.9602 do nbin=1,num_energybins e = ethresh-.125+nbin*0.25-mu energyprofile(nbin)=exp(-e*e/2/sigma/sigma) enddo * R3 distribution for PMT beta-gammas do rbin=1,num_r3bins r3 = (1.0*rbin)/num_r3bins*fidvol r3profile(rbin) = 1.631 + exp(-4.538 + 7.131*r3) enddo normalization=0. do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins component(nbin,rbin,cbin,zbin)= + energyprofile(nbin)*r3profile(rbin)*dnprofile(zbin) normalization = normalization + +component(nbin,rbin,cbin,zbin) enddo enddo enddo enddo do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins if (normalization .gt. 0) + pdf(nbin,rbin,cbin,zbin) = pdf(nbin,rbin,cbin,zbin) + +component(nbin,rbin,cbin,zbin) + *pmt_bgs/normalization enddo enddo enddo enddo * Now calculate the AV beta-gamma component * Energy spectrum of beta-gammas is tail of Gaussian, according to Vadim dnprofile(1) = settings(5,1) dnprofile(2) = settings(5,2) mu = 3.441 sigma = 0.4617 do nbin=1,num_energybins e = ethresh-.125+nbin*0.25-mu energyprofile(nbin)=exp(-e*e/2/sigma/sigma) enddo * R3 distribution for AV beta-gammas do rbin=1,num_r3bins r3 = (1.0*rbin)/num_r3bins*fidvol r3profile(rbin) = exp(-((1.056-r3)**2)/2/.1267/.1267) enddo normalization=0. do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins component(nbin,rbin,cbin,zbin)= + energyprofile(nbin)*r3profile(rbin)*dnprofile(zbin) normalization = normalization + +component(nbin,rbin,cbin,zbin) enddo enddo enddo enddo do nbin=1,num_energybins do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins if (normalization .gt. 0) + pdf(nbin,rbin,cbin,zbin) = pdf(nbin,rbin,cbin,zbin) + +component(nbin,rbin,cbin,zbin) + *av_bgs/normalization enddo enddo enddo enddo * Write out energy PDF of low E background do nbin=1,num_energybins normalization=0 do rbin=1,num_r3bins do cbin=1,num_cstsunbins do zbin=1,used_dnbins normalization=normalization+pdf(nbin,rbin,cbin,zbin) enddo enddo enddo if(normalization.ne.0.)write(*,*)'Low E',nbin,normalization *** alainb debug call hfill(4001,float(nbin)-0.5,0.,real(normalization)) *** alainb debug enddo * Write out R**3 PDF of low E background do rbin=1,num_r3bins normalization=0 do nbin=1,num_energybins do cbin=1,num_cstsunbins do zbin=1,used_dnbins normalization=normalization+pdf(nbin,rbin,cbin,zbin) enddo enddo enddo if(normalization.ne.0.)write(*,*)'R**3',rbin,normalization *** alainb debug call hfill(4002,float(rbin)-0.5,0.,real(normalization)) *** alainb debug enddo *** alainb debug * Write out costhetasun PDF of low E background do cbin=1,num_cstsunbins normalization=0 do nbin=1,num_energybins do rbin=1,num_r3bins do zbin=1,used_dnbins normalization=normalization+pdf(nbin,rbin,cbin,zbin) enddo enddo enddo if(normalization.ne.0.)write(*,*)'CTHSUN',cbin,normalization call hfill(4003,float(cbin)-0.5,0.,real(normalization)) enddo call hrput(0,'alainb.hist','n') *** alainb debug return end function azimuth (drxsun, drysun) implicit none real drxsun, drysun, azimuth real b, c real pi /3.1415926536/ b = drxsun/(drxsun**2)**0.5 c = (drxsun**2+drysun**2)**0.5 azimuth = b*acos(drysun/c)+pi*(1-b) - .86533 if (azimuth.le.0) then azimuth = azimuth+2*pi endif return end