REAL FUNCTION myheepcheck(a) IMPLICIT NONE INCLUDE '?' LOGICAL FIRST DATA FIRST /.TRUE./ SAVE FIRST REAL xcointime REAL p_e, th_e, p_p, th_p,newE, phi_e, phi_p REAL Hdist, Aldist, myhstheta, mysstheta, myhsphi, myssphi REAL hs_spec_angle, ss_spec_angle REAL myEm, myw, Qsq * corsi means corrected singes, corsi_e calculates the what should be * seen in the proton arm using the electron and then subtracts PRENERGY * Both of these should be centered around zero REAL mycorsi_e, mycorsi_p REAL qx,qy,qz,myqabs,mypmx,mypmy,mypmz,mypm REAL myPmpar,myPmper,mypmoop REAL ekinrec,m_rec REAL BEAMENERGY,omeg REAL pbeta,ebeam,cosgamma REAL ethick,pthick,costheta,ethick_dummy, pthick_dummy REAL loss_13,loss_21,loss_22,loss_23 REAL loss_02,loss_11,loss_12,loss_01 REAL ELENERGY,PRENERGY REAL m_elec, m_prot, m_amu, tar_a, tar_amin1, tar_z, tar_dens REAL e_beam, hs_spec_angle_deg, ss_spec_angle_deg REAL tar_thick, tar_radius REAL dE, dp_ss, dth_ss, dphi_ss, dp_hs, dth_hs, dphi_hs, set_sos_e REAL rangle, rqx, rqy, rqz, rmypmx, rmypmy, rmypmz REAL newhsxptar, newssxptar REAL rcW,rcEm,rcPMpar,rcPmper,rcPmoop,rcecorsi,rcpcorsi,rcphi REAL newssyptar, newhsyptar, phi_check logical good_electron logical good_proton logical coint_e_p REAL a(22) e_beam = a(1) set_sos_e = a(2) hs_spec_angle_deg = a(3) ss_spec_angle_deg = a(4) xcointime = a(5) dE= a(6) newssyptar=ssyptar-a(7)/1000.0 dp_ss= a(8) newhsyptar=hsyptar+a(9)/1000.0 dp_hs= a(10) dphi_ss= a(11) dphi_hs= a(12) newhsxptar= hsxptar- a(13) newssxptar= ssxptar- a(14) rcW= a(15) rcEm= a(16) rcPmpar= a(17) rcPmper= a(18) rcPmoop= a(19) rceCorsi= a(20) rcpCorsi= a(21) rcphi= a(22) m_elec=0.000511 m_prot=0.93827203 m_amu =0.93149404 tar_a = 1. tar_amin1 = 0. tar_z = 1. tar_dens=0.0723 tar_thick=0.288 tar_radius=2.0 * calculate lab frame angles from spectrometer gradients hs_spec_angle= hs_spec_angle_deg*3.14159265358/180.0 cosgamma=1.0/sqrt(1.0+newhsxptar**2+newhsyptar**2) costheta=cosgamma*( cos(hs_spec_angle) + $ newhsyptar*sin(hs_spec_angle) ) myhstheta=acos(costheta) myhsphi=atan2(newhsxptar, $ (sin(hs_spec_angle)-newhsyptar*cos(hs_spec_angle)) ) ss_spec_angle= ss_spec_angle_deg*3.14159265358/180.0 cosgamma=1.0/sqrt(1.0+newssxptar**2+newssyptar**2) costheta=cosgamma*( cos(ss_spec_angle) - $ newssyptar*sin(ss_spec_angle) ) mysstheta=acos(costheta) myssphi=atan2(newssxptar, $ (-sin(ss_spec_angle)-newssyptar*cos(ss_spec_angle)) ) if (myssphi<-1.5) then !remove the branch cut myssphi=myssphi+2.0*3.14159265358 endif * tweak momentum by an offest and set up variables with particle names if (set_sos_e.gt.0.5) then !electron in the sos * print *, 'electron in the SOS' p_e = ssp*(1-dp_ss/1000.0) th_e = mysstheta phi_e = myssphi p_p = hsp*(1-dp_hs/1000.0) th_p = myhstheta phi_p = myhsphi * phi_check=p_e*newssxptar+p_p*newhsxptar phi_check=newssxptar+(p_p/p_e)*newhsxptar else * print *, 'electron in HMS' p_e = hsp*(1-dp_hs/1000.0) th_e = myhstheta phi_e = myhsphi p_p = ssp*(1-dp_ss/1000.0) th_p = mysstheta phi_p = myssphi * phi_check=p_p*newssxptar+p_e*newhsxptar phi_check=(p_p/p_e)*newssxptar+newhsxptar endif newE=e_beam*(1+dE/1000.0) PRENERGY = sqrt(p_p**2 + m_prot**2) ELENERGY = p_e BEAMENERGY = newE *these components are in the lab frame qx = -p_e*sin(th_e)*cos(phi_e) ! q momentum is opposite to electron qy = -p_e*sin(th_e)*sin(phi_e) qz = BEAMENERGY - p_e*cos(th_e) mypmx = qx - p_p*sin(th_p)*cos(phi_p) mypmy = qy - p_p*sin(th_p)*sin(phi_p) mypmz = qz - p_p*cos(th_p) *Define new basis, the lab frame rotated about z, *so that q_vec is in the rotated xz plane rangle=phi_e rqx= qx*cos(rangle)+qy*sin(rangle) rqy= -qx*sin(rangle)+qy*cos(rangle) rqz = qz rmypmx= mypmx*cos(rangle)+mypmy*sin(rangle) rmypmy= -mypmx*sin(rangle)+mypmy*cos(rangle) rmypmz = mypmz myqabs= sqrt(rqx**2+rqy**2+rqz**2) myPm = sqrt(rmypmx**2 + rmypmy**2 + rmypmz**2) if(myqabs.ne.0.0)then myPmpar = (rmypmx*rqx+rmypmy*rqy+rmypmz*rqz)/myqabs myPmper = (-rqz*rmypmx+rqx*rmypmz)/myqabs myPmoop = rmypmy else myPmpar = 0.0 myPmper=0.0 myPmoop = 0.0 endif omeg = BEAMENERGY-p_e if(tar_amin1.ge.999. .or.tar_amin1.lt.0.001)then !Not used ekinrec = 0 else ekinrec = sqrt(myPm**2 + m_rec**2)-m_rec endif myEm= omeg +m_prot-PRENERGY - ekinrec Qsq=4.0*BEAMENERGY*ELENERGY*(sin(th_e/2.0))**2 myW = sqrt(m_prot**2+2.*m_prot*omeg-Qsq) mycorsi_e=p_e-BEAMENERGY/ + (1.+(BEAMENERGY/m_prot)*(1-cos(th_e))) mycorsi_p=p_p-2.*m_prot*BEAMENERGY*cos(th_p)/(BEAMENERGY+m_prot)/ + (1.-(BEAMENERGY*cos(th_p)/(BEAMENERGY+m_prot))**2) ******************************************************************* * Fill histograms using cuts if (set_sos_e.gt.0.5) then c ELECTRON IS IN THE SOS good_electron = .false. if ( ^ scer_npe .gt. 1.0 ^ .and.abs(ssyptar).lt.0.07 !next three are acceptance cuts ^ .and.abs(ssxptar).lt.0.02 ^ .and.abs(ssdelta-2.5).lt.12.5 ^ ) ^ good_electron = .true. c proton is in the HMS, check if it is an elastic proton good_proton = .false. if ( ^ hcer_npe.lt.1.0 ^ .and.abs(hsyptar).lt.0.04 !next three are acceptance cuts ^ .and.abs(hsxptar).lt.0.07 ^ .and.abs(hsdelta).lt.8.0 ^ .and.abs(hsztar-0.6).lt.1.2 !try to cut the target walls ^ ) ^ good_proton = .true. c now check for coincidence coint_e_p = .false. if( ^ good_electron ^ .and.good_proton ^ .and.abs(cointime-xcointime).lt.0.5 ^ ) ^ coint_e_p=.true. else c ELECTRON IS IN THE HMS good_electron = .false. if ( ^ hcer_npe .gt. 1.0 ^ .and.abs(hsyptar).lt.0.04 !next three are acceptance cuts ^ .and.abs(hsxptar).lt.0.07 ^ .and.abs(hsdelta).lt.8.0 ^ .and.abs(hsztar-0.6).lt.1.2 !try to cut the target walls ^ ) ^ good_electron = .true. c proton is in the SOS, check if it is an elastic proton good_proton = .false. if ( ^ scer_npe.lt.1.0 ^ .and.abs(ssyptar).lt.0.07 !next three are acceptance cuts ^ .and.abs(ssxptar).lt.0.02 ^ .and.abs(ssdelta-2.5).lt.12.5 ^ ) ^ good_proton = .true. c now check for coincidence coint_e_p = .false. if( ^ good_electron ^ .and.good_proton ^ .and.abs(cointime-xcointime).lt.0.5 ^ ) ^ coint_e_p=.true. endif * Create histograms IF(FIRST) THEN FIRST=.FALSE. PRINT *, 'BOOKING HISTOS...' CALL HBOOK1(3,'myW (MeV)-mprot',100, $ -200.0,200.0,0) CALL HBOOK1(4,'myEm (MeV)',100,-30.0,30.0,0) CALL HBOOK1(5,'myPmpar (MeV)',100,-100.0,100.0,0) CALL HBOOK1(6,'myPmper (MeV)',100,-100.0,100.0,0) CALL HBOOK1(7,'myPmoop (MeV)',100,-100.0,100.0,0) CALL HBOOK1(8,'myCorsi (e) (MeV)',100,-100.0,100.0,0) CALL HBOOK1(9,'myCorsi (p) (MeV)',100,-100.0,100.0,0) c CALL HBOOK1(51,'SOS phi-HMS phi-180 (deg)',100,-10.0,10.0,0) c CALL HBOOK1(51,'ssxptar*ssp+hsxptar*hsp',100,-0.1,0.1,0) CALL HBOOK1(51,'(e xptar)+(Pp/Pe)*(xptar p) (mrad)' $ ,100,-100.,100.,0) CALL HBARX(3) CALL HBARX(4) CALL HBARX(5) CALL HBARX(6) CALL HBARX(7) CALL HBARX(8) CALL HBARX(9) CALL HBARX(51) c$$$ CALL HBOOK2(51, 'bonus plot', 100, 1.5, c$$$ ^ 1.8, 100, 0.98, 1.02, 0.0) ENDIF * Fill histograms if(coint_e_p)then CALL HFILL(3, (myW-m_prot)*1000.0-rcW ,0.,1.) CALL HFILL(4, (myEm)*1000.0-rcEm ,0.,1.) CALL HFILL(5, (myPmpar)*1000.0-rcPmpar ,0.,1.) CALL HFILL(6, (myPmper)*1000.0-rcPmper ,0.,1.) CALL HFILL(7, (myPmoop)*1000.0-rcPmoop ,0.,1.) CALL HFILL(8, (myCorsi_e)*1000.0-rceCorsi ,0.,1.) CALL HFILL(9, (myCorsi_p)*1000.0-rcpCorsi ,0.,1.) * CALL HFILL(51, phi_check*1000.0 ,0.,1.) CALL HFILL(51, phi_check*1000.0 ,0.,1.) c CALL HFILL(51, 57.296*(myssphi-myhsphi)-180.0-rcphi ,0.,1.) c$$$ CALL hf2(51, mypmper,pmper, 1.0) endif return END