************************************************************************* * This version of momanhough performs tracker CALIBRATIONS: * - eta functions * - charge correlation * - ADC-to-MIP conversion * ------------------------------------------------------------------- * It works as it follow: * - it reads cluster information (read a list of LEVEL1 ntuples, * output of reduction.f) * - it performs track identification, select good tracks (6 hit planes) * and perform first-order fit * - it calculate calibration variables * - it creates a LEVEL2 ntuple with selected events (in order to check) * ************************************************************************* program momanhough include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' include '../common/common_level2debug.f' include '../common/common_mech.f' include '../common/common_xyzPAM.f' include '../common/common_mini_2.f' include '../common/calib.f' include '../common/level1.f' include '../common/level2.f' include '../common/common_calibration.f' * flag to set debug mode logical DEBUG common/dbg/DEBUG logical DONTDO data DONTDO/.true./ c$$$* flag to chose PFA c$$$ character*10 PFA c$$$ common/FINALPFA/PFA c parameter (inf=1.e8) !just a huge number... * external functions external npl external acoordsi,coordsi,nld,coord,dcoord c------------------------------------------------------------------------ c c local variables c c------------------------------------------------------------------------ character*24 processing_date parameter (nfile_max=500) parameter (lun_data_level1=71) !data file id number parameter (lun_data_level2=72) !data file id number parameter (lun_data_calib=74) !data file id number c parameter (lun_file_out=75) !data file id number character*74 data_file !data file name character*74 data_dir ! character*74 out_dir ! character*74 data_file_calib character*74 data_file_level1 character*74 data_file_level2 character*10 input_string(nfile_max) character*74 file_out character*74 fname_param character*7 rzdir # ifndef TEST2003 parameter(ncalibmax=50) character*40 file_calib(ncalibmax) parameter(lun_calib_list=66) !calibration list file id integer which_calib_last # endif integer tottr integer minevent !first event to be analysed logical FIMAGE ! COMMON/QUEST/IQUEST(100) c !permette di ottenere ntuple funzionanti nonostante c ! il messaggio dei 64K di RZOUT...!??? # ifdef TEST2003 parameter (nang=3) !number of angular bins parameter (angstep=4.) parameter (angmin=-6.) print*,' ' print*,'Set angular binning TEST 2003' print*,' ' nangbin=nang do ia=1,nangbin angL(ia)=angmin+angstep*(ia-1) angR(ia)=angmin+angstep*ia print*,ia,' - ',angL(ia),angR(ia) enddo print*,' ' # else c$$$ parameter (nang=21) !number of angular bins c$$$ parameter (angstep=2.) c$$$ parameter (angmin=-21.) c$$$ print*,' ' c$$$ print*,'Set angular binning Rome-2005' c$$$ print*,' ' c$$$ c$$$ nangbin=nang c$$$ do ia=1,nangbin c$$$ angL(ia)=angmin+angstep*(ia-1) c$$$ angR(ia)=angmin+angstep*ia c$$$ print*,ia,' - ',angL(ia),angR(ia) c$$$ enddo c$$$ print*,' ' parameter (nang=17) !number of angular bins * NB angL and angR are created with length 21 data angL $ /-21.,-17.,-13.,-11.,-9.,-7.,-5.,-3.,-1. $ ,1.,3.,5.,7.,9.,11.,13.,17.,0.,0.,0.,0./ data angR $ /-17.,-13.,-11.,-9.,-7.,-5.,-3.,-1. $ ,1.,3.,5.,7.,9.,11.,13.,17.,21.,0.,0.,0.,0./ print*,' ' print*,'Set angular binning Rome-2005' print*,' ' nangbin=nang do ia=1,nangbin print*,ia,' - ',angL(ia),angR(ia) enddo print*,' ' # endif c------------------------------------------------------------------------ c c HBOOK initialization c c------------------------------------------------------------------------ call HLIMIT(NWPAWC) c------------------------------------------------------------------------ c c reads input informations c c------------------------------------------------------------------------ call fdate(processing_date) write(*,101) $ processing_date 101 format(/ $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ $ ,'* *',/ $ ,'* CALIBRATION *',/ $ ,'* *',/ $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ $ ,a24,/ $ ) print*,'Number of files to be analysed:' read(*,*)nfile print*,nfile 111 format(a) print*,'Input data directory:' read(*,111)data_dir print*,data_dir print*,'Output data directory:' read(*,111)out_dir print*,out_dir print*,'List of file identifiers: (DATE_NUM)' 300 format(A10) do i=1,nfile read(*,300)input_string(i) print*,input_string(i) enddo minevent=1 print*,'Maximum number of events to be analized:' !20000 read(*,*) ntotev print*,ntotev print*,'DEBUG mode: (T, F)' 11 format(l1) read(*,11)DEBUG print*,DEBUG print*,'---------------------------------------------------' ****** INITIALIZATIONS ************************************* * 2) read z coordinates of the planes print*,' ' print*,'- read z coordinates of the planes' print*,' ' call mech_sensor !reads sensors centres coordinates do ip=1,nplanes fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm * gets planes mechanical z positions * (in mm) and sets them in micrometers enddo * 4) read magnetic field map print*,' ' print*,'- read magnetic field map' print*,' ' call read_B * 5) read allignment parameters print*,' ' print*,'- read aligment parameters' print*,' ' call readalignparam ************************************************************ ************************************************************ ************************************************************ ************************************************************ print*,' ' print*,'---------------------------------------------------' print*,' ' print*,'OPENING OUTPUT-HISTO FILES:' 1010 format('LADDER',i1) do ilad=1,3 write(rzdir,1010)ilad file_out= $ out_dir(1:LNBLNK(out_dir))//'trk-' $ //rzdir//'-histos.rz' print*,' -> '//file_out IQUEST(10)=65000 call HROPEN(lun_file_out+ilad,rzdir,file_out,'PQN',4096,istat) if(istat.ne.0) goto 16 call HCDIR('//'//rzdir,' ') call HCDIR('//PAWC','') call HMDIR(rzdir,'S') call book_eta_histos call book_charge_histos call book_check_histos call HCDIR('//PAWC','') enddo ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ ************************************************************ c------------------------------------------------------------------------ c c LEVEL 2 ntuple booking c c------------------------------------------------------------------------ 503 format(a,'DW__level2-calib.rz') write(data_file_level2,503) $ out_dir(1:LNBLNK(out_dir)) c $ ,data_file(1:LNBLNK(data_file)) print*,'' print*,'------------------------------------' print*,' Creating LEVEL2 rz file' print*,' ',data_file_level2 print*,'------------------------------------' C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C largest RZ file: IQUEST(10) records x LREC words x 4 byte C with the following settings: 65000 x 4096 x 4 = 1G C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IQUEST(10)=65000 c !permette di ottenere ntuple funzionanti nonostante c ! il messaggio dei 64K di RZOUT...!??? call HROPEN(lun_data_level2, $ 'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz if(istat.ne.0) goto 19 call book_level2 print*,' ' print*,'---------------------------------------------------' print*,' * * * START ANALYSIS * * *' print*,'---------------------------------------------------' c------------------------------------------------------------------------ c c loop on files c c------------------------------------------------------------------------ goodtr = 0 tottr = 0 do ifile = 1,nfile !files loop print*,' ' data_file = input_string(ifile) print*,'--> FILE ',ifile,': '//data_file(1:LNBLNK(data_file)) c------------------------------------------------------------------------ c c fills bad variables with online DSP bad strips from calib histograms c for bad strip exclusion c c------------------------------------------------------------------------ # ifdef TEST2003 c------------------------------------------------------------------------ c c opens calib file c c------------------------------------------------------------------------ 505 format(a,'output_',a,'_calib.rz') write(data_file_calib,505) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,' ' print*,'OPENING CALIB FILE...', data_file_calib IQUEST(10)=65000 call HROPEN $ (lun_data_calib,'IN',data_file_calib,'QP',4096,istat) if(istat.ne.0) goto 17 call HCDIR('//IN',' ') call HRIN(0,9999,0) !puts histograms in memory print*,' ' print*,' ' print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...' print*,' ' call fillpedsig call HREND('IN') close (lun_data_calib) # else print*,' ' print*,'OPENING CALIBRATION-LIST FILE:' 501 format(a,'DW_',a,'_calib.txt') write(data_file_calib,501) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,data_file_calib open(lun_calib_list, $ FILE=data_file_calib(1:LNBLNK(data_file_calib)) $ ,STATUS='UNKNOWN' $ ,IOSTAT=iostat $ ) 113 format(i5,' ',a25) do i=1,ncalibmax read(lun_calib_list,113,IOSTAT=iostat) n_cal_list $ ,file_calib(i) c file_calib(i)=file_calib(i)(1:LNBLNK(file_calib(i))) if(iostat.ne.0)then ncal=i-1 goto 2000 endif print*,n_cal_list,' - ' $ ,file_calib(i)(1:LNBLNK(file_calib(i))) enddo 2000 close(lun_calib_list) which_calib_last=0 # endif c------------------------------------------------------------------------ c c opens level1 file c c------------------------------------------------------------------------ # ifdef TEST2003 504 format(a,'output_',a,'_level1.rz') # else 504 format(a,'DW_',a,'_level1.rz') # endif write(data_file_level1,504) $ data_dir(1:LNBLNK(data_dir)) $ ,data_file(1:LNBLNK(data_file)) print*,'' print*,'OPENING LEVEL1 FILE:' print*,data_file_level1 IQUEST(10)=65000 c !permette di ottenere ntuple funzionanti nonostante c ! il messaggio dei 64K di RZOUT...!??? call HROPEN(lun_data_level1, $ 'LEVEL1',data_file_level1,'QP',4096,istat) !opens rz if(istat.ne.0) goto 19 call HRIN(ntp_level1,9999,20) call access_level1 c call HPRNTU(ntp_level1+20) call HNOENT(ntp_level1+20,iemax0) print*,'( ',iemax0,' events)' ************************************************************ ************************************************************ ************************************************************ * * start track analysis * ************************************************************ ************************************************************ ************************************************************ nselect = 0 maxevent=minevent+ntotev-1 do iev = minevent,MIN(iemax0,maxevent) !loop on events call HCDIR('//LEVEL1',' ') call HGNT(ntp_level1+20,iev,ierr) !reads an event if(ierr.ne.0) goto 21 tottr = tottr +1 *------------------------------------------------------ * LEVEL2 N-TUPLE INITIALIZATIONS call init_level2 if(.not.good1)then * goto 8800 !fill nt-uple and go to next event goto 100 !go to next event endif *------------------------------------------------------ # ifndef TEST2003 if(which_calib.ne.which_calib_last.and. $which_calib.ne.0)then data_file_calib= $ data_dir(1:LNBLNK(data_dir))// $ file_calib(which_calib) $ (1:LNBLNK(file_calib(which_calib))) c print*,data_file_calib print*,'' print*, $ '@ event ',nev2 $ ,' (CPU pkt N.',pkt_num1,')' print*,'--> ',data_file_calib IQUEST(10)=65000 call HROPEN(lun_data_calib, $ 'CALIB',data_file_calib,'QP',4096,istat) !opens if(istat.ne.0) goto 19 call HRIN(0,9999,0) call fillpedsig do iview=1,nviews call HDELET(id_hi_bad+iview) call HDELET(id_hi_ped+iview) call HDELET(id_hi_sig+iview) enddo call HREND('CALIB') close(lun_data_calib) which_calib_last=which_calib elseif(which_calib.eq.0)then nocalib=nocalib+1 good2=.false. goto 8800 !fill nt-uple and go to next event endif # endif *------------------------------------------------------ * cut on maximum number of clusters *------------------------------------------------------ * if(nclstr1.gt.nclstrmax_level2)then if(nclstr1.gt.20)then * goto 8800 !fill nt-uple and go to next event goto 100 !go to next event endif do i=1,nclstr1 cl_used(i)=0 !init mask of clusters associated to a track enddo if(DEBUG)then print*,'----------------------------------' print*,iev,' ** ',nev2 endif *------------------------------------------------------------------------------- * STEP 1 *------------------------------------------------------------------------------- * X-Y cluster association * * Clusters are associated to form COUPLES * Clusters not associated in any couple are called SINGLETS * * Track identification (Hough transform) and fitting is first done on couples. * Hence singlets are possibly added to the track. * * Variables assigned by the routine "cl_to_couples" are those in the * common blocks: * - common/clusters/cl_good * - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2 * - common/singlets/ncls,cls,cl_single * * NB - The routine version used for calibration is different from that used * for event processing. Different cluster selection is applied: * - "ngoodstr" good strip around MAXS required (ngoodstr=2 => 5 good strips) * - no charge correlation required *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- iflag=0 call cl_to_couples_nocharge(iflag) if(iflag.eq.1)then !bad event * goto 880 !fill ntp and go to next event goto 100 !go to next event endif *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PRE-SELEZIONE DI TRACCE PULITE *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * at least a couple per plane do ip=1,nplanes if(ncp_plane(ip).eq.0)goto 100 !next event enddo c if(good2.eq.0)goto 880!fill ntp and go to next event *----------------------------------------------------- *----------------------------------------------------- * HOUGH TRASFORM *----------------------------------------------------- *----------------------------------------------------- *------------------------------------------------------------------------------- * STEP 2 *------------------------------------------------------------------------------- * * Association of couples to form * - DOUBLETS in YZ view * - TRIPLETS in XZ view * * Variables assigned by the routine "cp_to_doubtrip" are those in the * common blocks: * - common/hough_param/ * $ alfayz1, !Y0 * $ alfayz2, !tg theta-yz * $ alfaxz1, !X0 * $ alfaxz2, !tg theta-xz * $ alfaxz3 !1/r * - common/doublets/ndblt,cpyz1,cpyz2 * - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3 *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- iflag=0 call cp_to_doubtrip(iflag) if(iflag.eq.1)then !bad event * goto 880 !fill ntp and go to next event goto 100 !go to next event endif *------------------------------------------------------------------------------- * STEP 3 *------------------------------------------------------------------------------- * * Classification of doublets and triplets to form CLOUDS, * according to distance in parameter space. * * cloud = cluster of points (doublets/triplets) in parameter space * * * * Variables assigned by the routine "doub_to_YZcloud" are those in the * common blocks: * - common/clouds_yz/ * $ nclouds_yz * $ ,alfayz1_av,alfayz2_av * $ ,ptcloud_yz,db_cloud,cpcloud_yz * * Variables assigned by the routine "trip_to_XZcloud" are those in the * common blocks: * common/clouds_xz/ * $ nclouds_xz * $ ,alfaxz1_av,alfaxz2_av,alfaxz3_av * $ ,ptcloud_xz,tr_cloud,cpcloud_xz *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- iflag=0 call doub_to_YZcloud(iflag) if(iflag.eq.1)then !bad event * goto 880 !fill ntp and go to next event goto 100 !go to next event endif iflag=0 call trip_to_XZcloud(iflag) if(iflag.eq.1)then !bad event * goto 880 !fill ntp and go to next event goto 100 !go to next event endif *------------------------------------------------------------------------------- * STEP 4 (ITERATED until any other physical track isn't found) *------------------------------------------------------------------------------- * * YZ and XZ clouds are combined in order to obtain the initial guess * of the candidate-track parameters. * A minimum number of matching couples between YZ and XZ clouds is required. * * A TRACK CANDIDATE is defined by * - the couples resulting from the INTERSECTION of the two clouds, and * - the associated track parameters (evaluated by performing a zero-order * track fitting) * * The NTRACKS candidate-track parameters are stored in common block: * * - common/track_candidates/NTRACKS,AL_STORE * $ ,XV_STORE,YV_STORE,ZV_STORE * $ ,XM_STORE,YM_STORE,ZM_STORE * $ ,RESX_STORE,RESY_STORE * $ ,AXV_STORE,AYV_STORE * $ ,XGOOD_STORE,YGOOD_STORE * $ ,CP_STORE,RCHI2_STORE * *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- ntrk=0 !counter of identified physical tracks 11111 continue !<<<<<<< come here when performing a new search iflag=0 call clouds_to_ctrack(iflag) if(iflag.eq.1)then !no candidate tracks found * goto 880 !fill ntp and go to next event goto 100 !go to next event endif *------------------------------------------------------ * cut on maximum number of track candidates * (to have clean events) * NB: must be at least 2 to do not discard tracks * contained in a sensor column on Y view * * cut replaced with cut on the # cluster out of the track *------------------------------------------------------ * if(ntracks.gt.2)then * goto 8800 !fill nt-uple and go to next event * goto 100 !go to next event * endif * if(DEBUG)print*,'*GOOD!*' FIMAGE=.false. !processing best track (not track image) ibest=0 !best track among candidates iimage=0 !track image * ----------------------------------------------------- * --------------- select the best track --------------- * selection cuts: * - NX = NY = 6 * ----------------------------------------------------- * rchi2best=100000. do i=1,ntracks npx=0 npy=0 do ip=1,nplanes npx=npx+XGOOD_STORE(ip,i) npy=npy+YGOOD_STORE(ip,i) enddo * print*, * $ 'Candidate ',i,' nx,ny ',npx,npy,' chi2 ',RCHI2_STORE(i) if(npx.eq.nplanes.and.npy.eq.nplanes)then if(RCHI2_STORE(i).lt.rchi2best.and. $ RCHI2_STORE(i).gt.0)then ibest=i rchi2best=RCHI2_STORE(i) endif endif enddo * if(ibest.eq.0)goto 880 !>> no good candidates if(ibest.eq.0)goto 100 !>> no good candidates * print*,'BEST TRACK ',ibest,' nx,ny',npx,npy * print*,'>>>> ',iev * ----------------------------------------------------- * ------------- store the best track only ------------- * ----------------------------------------------------- ntr=1 ntrk=ntr image(ntr)=0 good2=.true. chi2_nt(ntr) = RCHI2_STORE(ibest) c print*,chi2_nt(ntr) * **************************************************** * normalize the angles * **************************************************** pig=acos(-1.) phi = AL_STORE(4,ibest) sinth = AL_STORE(3,ibest) if(sinth.lt.0)then sinth = -sinth phi = phi + pig endif npig = aint(phi/(2*pig)) phi = phi - npig*2*pig if(phi.lt.0) $ phi = phi + 2*pig AL_STORE(4,ibest) = phi AL_STORE(3,ibest) = sinth * **************************************************** * ... then assign the other variables * **************************************************** do i=1,5 al_nt(i,ntr) = AL_STORE(i,ibest) c print*,i,'al_nt(i,ntr)',al_nt(i,ntr) do j=1,5 coval(i,j,ntr) = 0 ! TEMPORANEO!!!!!!!!!!!! enddo enddo do ip=1,nplanes ! loop on planes xgood_nt(ip,ntr) = XGOOD_STORE(ip,ibest) ygood_nt(ip,ntr) = YGOOD_STORE(ip,ibest) xm_nt(ip,ntr) = XM_STORE(ip,ibest) ym_nt(ip,ntr) = YM_STORE(ip,ibest) zm_nt(ip,ntr) = ZM_STORE(ip,ibest) RESX_nt(IP,ntr) = RESX_STORE(ip,ibest) RESY_nt(IP,ntr) = RESY_STORE(ip,ibest) xv_nt(ip,ntr) = XV_STORE(ip,ibest) yv_nt(ip,ntr) = YV_STORE(ip,ibest) zv_nt(ip,ntr) = ZV_STORE(ip,ibest) axv_nt(ip,ntr) = AXV_STORE(ip,ibest) ayv_nt(ip,ntr) = AYV_STORE(ip,ibest) c print*,ip,' xm_nt(ip,ntr)', xm_nt(ip,ntr) c print*,ip,' ym_nt(ip,ntr)', ym_nt(ip,ntr) c print*,ip,' zm_nt(ip,ntr)', zm_nt(ip,ntr) c print*,ip,' RESX_nt(ip,ntr)', RESX_nt(ip,ntr) c print*,ip,' RESY_nt(ip,ntr)', RESY_nt(ip,ntr) c print*,ip,' xv_nt(ip,ntr)',xv_nt(ip,ntr) c print*,ip,' yv_nt(ip,ntr)',yv_nt(ip,ntr) c print*,ip,' zv_nt(ip,ntr)',zv_nt(ip,ntr) c print*,ip,' axv_nt(ip,ntr)',axv_nt(ip,ntr) c print*,ip,' ayv_nt(ip,ntr)',ayv_nt(ip,ntr) id = CP_STORE(ip,ibest) icp = icp_cp(id) iclx = clx(nplanes-ip+1,icp) icly = cly(nplanes-ip+1,icp) which_cl_x(ip) = iclx which_cl_y(ip) = icly dedx_x(ip,ntr) = dedx(iclx) dedx_y(ip,ntr) = dedx(icly) enddo c call CalcBdL(100,xxxx,IFAIL) c if(ifps(xxxx).eq.1)BdL(ntr) = xxxx BdL(ntr) = 0 * --- then remove selected clusters (ibest+iimage) from clouds ---- call clean_XYclouds(ibest) if(iflag.eq.1)then !bad event * goto 880 !fill ntp and go to next event goto 100 !go to next event endif * >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<< * + + + + + + + + + + + + + + + + + * / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / * + + + + + + + + + + + + + + + + + 880 continue * ********************************************************** * stores info about clusters not associated with any track * ********************************************************** c call fill_level2_siglets if(DEBUG)then print*,'' print*,'DONE!' print*,'' print*,'* summary *' print*,'tracks ',ntrk print*,'cl used ',(cl_used(i),i=1,nclstr1) c$$$ print*,'cl unused (x-y)' c$$$ do ip=1,nplanes c$$$ print*,ip,' << ',nclsx(ip),nclsy(ip) c$$$ enddo print*,'' print*,'' endif 8800 continue * ----------------------------------------------------- * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * calibrations: * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * ----------------------------------------------------- * ================================== * ////////////////////////////////// * Track selection * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ * ================================== * - cut on number of planes * (fatto a monte... non mi ricordo perche') c$$$ npx=0 c$$$ npy=0 c$$$ do ip=1,nplanes c$$$ npx=npx+XGOOD_STORE(ip,i) c$$$ npy=npy+YGOOD_STORE(ip,i) c$$$ enddo c$$$ if(.not.(npx.ge.5..and.npy.ge.5))goto 100 * - chi^2 selection rig = abs(1./al_nt(5,1)) c$$$ chicut = 5. !100. c$$$ chicut = chicut * ( 1.48 + 4.64/rig + 18.35/rig**2 ) c$$$ if(chi2_nt(1).gt.chicut)goto 100 !next event chicut = 20 chicut = chicut *( c $ 3.9633*exp(-0.54485E-01*rig) c $ +0.99231 c $ +10.835*exp(-0.67155*rig) ) $ 3.7032*exp(-0.41911E-01*rig) $ +0.85355 $ +10.905*exp(-0.61893*rig) ) if(chi2_nt(1).gt.chicut.or.chi2_nt(1).lt.0)goto 100 !next event * - cut on number of cluster out of the track (noise) * (no constraints on the last plane becouse of the * back-scattering from calorimeter) c$$$ do ip=1,5 !6 c$$$ nclstot=nclsx(ip)+nclsy(ip) c$$$ if(nclstot.ge.1)goto 100 !next event c$$$ enddo nclstot = 0 do isx=1,nclsx c if(planex(isx).ne.6)nclstot = nclstot +1 nclstot = nclstot +1 enddo do isy=1,nclsy c if(planey(isy).ne.6)nclstot = nclstot +1 nclstot = nclstot +1 enddo if(nclstot.ne.0)print*,'iev ',iev,' ncls ',nclstot c if(nclstot.ge.1)goto 100 * ================================== * ////////////////////////////////// * Track selection (end) * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ * ================================== call HCDIR('//LEVEL2',' ') goodtr = goodtr +1 nselect = nselect +1 c print*,goodtr,' ****',iev call HFNT(ntp_level2) !fill LEVEL2 nt-uple call HCDIR('//PAWC','') * ----------------------------------------------------- * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * calibrations (sum over events) * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * ----------------------------------------------------- if(goodtr.lt.ncorrmax)then do ip=1,nplanes !1-->6 iviewx = (nplanes - ip)*2+2 !12-->2 ilx=nld(maxs(which_cl_x(ip)),iviewx) iviewy = (nplanes - ip)*2+1 !11-->1 ily=nld(maxs(which_cl_y(ip)),iviewy) if(ilx.ne.ily)print*,'LADDERS DO NOT MATCH!!!' chargexy(ip,ilx,goodtr,1)=dedx(which_cl_x(ip)) chargexy(ip,ily,goodtr,2)=dedx(which_cl_y(ip)) enddo endif * compute average amgle in each bin do i=1,nviews do ia=1,nangbin if(navangbin(ia,i).ne.0) $ avangbin(ia,i)= $ avangbin(ia,i)/navangbin(ia,i) enddo enddo * apply a cut on the measured rigidity: if(rig.gt.rigmax)then call fill_eta_histos call fill_check_histos endif * ----------------------------------------------------- * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * calibrations (end) * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * ----------------------------------------------------- 100 continue enddo !end loop on events print*,'' print *,'- Stored ',nselect,' entries in LEVEL2 nt-uple' call HCDIR('//LEVEL1',' ') call HREND('level1') close(lun_data_level1) call HCDIR('//PAWC',' ') call HDELET(ntp_level1+20) enddo !end loop on files print*,'' call HCDIR('//LEVEL2',' ') call RZPURG(-1) call HROUT(ntp_level2,ICYCLE,'T') call HREND('level2') close(lun_data_level2) call HCDIR('//PAWC',' ') call HDELET(ntp_level2) * =================================================== * in the following the charge correlation parameters * and the eta correction functions are evaluated * and stored in txt files. * However I chosed to do it "off-line", by means of * a kumac, so that I can group more files together. * =================================================== * ----------------------------------------------------- * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * calibrations * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * ----------------------------------------------------- do il=1,3 write(rzdir,1010)il ************************************************** ************************************************** * CHARGE CORRELATION CALIBRATION ************************************************** ************************************************** if(DONTDO)goto 8881 call charge_calibration(il) 107 format(/ $ ,/ $ ,/ $ ' >>> CHARGE-CORRELATION CALIBRATION <<<',/ $ ,/ $ 'Save charge correlation parameters --->',/ $ ,a,/ $ ,/ $ ,'Parameterization ------> Sy = K * Sx + C',/ $ ,'Plane K EK C EC',/ $ ,'-----------------------------------------') c write(lun_file_log,107)file_out(1:LNBLNK(file_out)) file_out= $ out_dir(1:LNBLNK(out_dir))//'trk-' $ //rzdir//'-charge.dat' write(*,107)file_out(1:LNBLNK(file_out)) open(55, $ FILE=file_out(1:LNBLNK(file_out)), c $ STATUS='NEW', $ FORM='FORMATTED') 108 format(i5,4f9.5) do ip=1,nplanes c write(lun_file_log,108) !log file c $ ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) write(*,108) !screen $ ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) write(55,*) !data $ ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) enddo close(55) 8881 call fill_charge_histos(il) ************************************************** ************************************************** * ETA CALIBRATION ************************************************** ************************************************** file_out= $ out_dir(1:LNBLNK(out_dir))//'trk-' $ //rzdir 109 format('Creating files:',/ $ ,' ',a,'-eta2-binning.dat',/ $ ) write(*,109)file_out(1:LNBLNK(file_out)) open(lun_file_eta, $ FILE=file_out(1:LNBLNK(file_out))//'-eta2-binning.dat' $ ,FORM='FORMATTED' $ ) do ia=1,nangbin write(lun_file_eta,*)ia,angL(ia),angR(ia) enddo close(lun_file_eta) if(DONTDO)goto 8882 do ia=1,nang 201 format('-eta2-bin',i1,'.dat') 202 format('-eta2-bin',i2,'.dat') if(ia.lt.10)write(fname_param,201)ia if(ia.ge.10)write(fname_param,202)ia open(lun_file_eta, $ FILE=file_out(1:LNBLNK(file_out)) $ //fname_param(1:LNBLNK(fname_param)) $ ,FORM='FORMATTED' $ ) call eta_calibration(il,ia) * ********************** write ************************** 203 format(13f10.4) do ie=1,nhbin2+1 write(lun_file_eta,203) $ eta2abs(ie),(eta2corr(ie,i),i=1,nviews) enddo close(lun_file_eta) * ********************** write ************************** enddo !end loop on angular bins 8882 continue enddo !end loop on LADDERs * ----------------------------------------------------- * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * calibrations * +++++++++++++++++++++++++++++++++++++++++++++++++++++ * ----------------------------------------------------- 8888 continue c------------------------------------------------------------------------ c c no error exit c c------------------------------------------------------------------------ goto 9000 !happy ending c------------------------------------------------------------------------ c c data file opening error c c------------------------------------------------------------------------ 19 continue print*,' ' print*,'ERROR OPENING DATA FILE: ',data_file print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c level1 ntuple event reading error c c------------------------------------------------------------------------ 21 continue print*,' ' print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT $ : ',iev print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c calib file opening error c c------------------------------------------------------------------------ 17 continue print*,' ' print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib print*,' ' print*,' ' goto 9000 !the end 16 continue print*,' ' print*,'preanalysis: ERROR OPENING OUTPUT FILE: ',file_out print*,' ' print*,' ' goto 9000 !the end c------------------------------------------------------------------------ c c closes files and exits c c------------------------------------------------------------------------ 9000 continue c$$$ if(DEBUG)call HPRNTU(ntp_level2+1) c$$$ call HPRNTU(ntp_level2) c$$$ print*,'' c$$$ call HCDIR('//LEVEL2',' ') c$$$ call HROUT(ntp_level2,ICYCLE,'T') c$$$ print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )' c$$$ if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T') c$$$ if(DEBUG)print *,'- Stored DEBUG nt-uple ' c$$$ call HREND('level2') c$$$ close(lun_data_level2) c$$$ call HCDIR('//LEVEL1',' ') c$$$ call HREND('level1') c$$$ close(lun_data_level1) write(*,105) $ tottr $ ,goodtr 105 format(/ $ ,/ $ ,'|------------------------------------------------|',/ $ ,'| Number of processed events |',/ $ ,'| Total ',i8,' |'/ $ ,'| Good tracks ',i8,' |'/ $ ,'|------------------------------------------------|',/ $ ) do ilad=1,3 write(rzdir,1010)ilad call HCDIR('//PAWC/'//rzdir,' ') call HCDIR('//'//RZDIR,' ') call HROUT(id_hi_chi,ICYCLE,'T') call HROUT(id_hi_chi+100,ICYCLE,'T') call HROUT(id_hi_chi+200,ICYCLE,'T') call HROUT(id_hi_rig,ICYCLE,'T') call HROUT(id_hi_charge,ICYCLE,'T') do i=1,nviews call HROUT(id_hi_charge+i*100,ICYCLE,'T') call HROUT(id_hi_charge+i*100+1,ICYCLE,'T') call HROUT(id_hi_charge+i*100+2,ICYCLE,'T') call HROUT(80000+id_hi_charge+i*100+2,ICYCLE,'T') call HROUT(id_hi_mult+i*100,ICYCLE,'T') call HROUT(80000+id_hi_mult+i*100,ICYCLE,'T') do ia=1,nang call HROUT(id_hi_eta2+i*100+ia,ICYCLE,'T') call HROUT(id_hi_eta3+i*100+ia,ICYCLE,'T') call HROUT(id_hi_eta4+i*100+ia,ICYCLE,'T') c call HROUT(id_hi_charge+i*100+ia,ICYCLE,'T') call HROUT(id_hi_mult+i*100+ia,ICYCLE,'T') call HROUT(id_hi_mult+i*100+ia+50000,ICYCLE,'T') enddo c call HROUT(id_hi_chi+100+ia,ICYCLE,'T') c call HROUT(id_hi_chi+200+ia,ICYCLE,'T') enddo do ip=1,nplanes call HROUT(id_hi_angle+100+ip,ICYCLE,'T') call HROUT(id_hi_angle+200+ip,ICYCLE,'T') call HROUT(id_hi_angle_stat+100+ip,ICYCLE,'T') call HROUT(id_hi_angle_stat+200+ip,ICYCLE,'T') call HROUT(id_hi_chargeco+ip*100+4,ICYCLE,'T') call HROUT(id_hi_chargeco+ip*100,ICYCLE,'T') do is=1,2 call HROUT(id_hi_beam+ip*100+is+200000,ICYCLE,'T') call HROUT(id_hi_beam+ip*100+is+100000,ICYCLE,'T') enddo call HROUT(id_hi_bx+ip,ICYCLE,'T') call HROUT(id_hi_by+ip,ICYCLE,'T') call HROUT(id_hi_bz+ip,ICYCLE,'T') enddo call HREND(RZDIR) close(lun_file_out+ilad) enddo stop end # include "momanhough-subroutines.F" ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** ******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** ******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** ******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** ******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** ******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** ******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* *************************************************************** * * * * * *************************************************************** subroutine book_eta_histos include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' character*64 title !histos title 502 format('Eta-2 (view ',i2,', angle bin ',i2,' )') 503 format('Eta-3 (view ',i2,', angle bin ',i2,' )') 504 format('Eta-4 (view ',i2,', angle bin ',i2,' )') do i=1,nviews do ia=1,nangbin write(title,502)i,ia call HBOOK1(id_hi_eta2+i*100+ia $ ,title c $ ,nhbin2,0.,1.,0.) $ ,nhbin2,-0.5,0.5,0.) write(title,503)i,ia call HBOOK1(id_hi_eta3+i*100+ia $ ,title c $ ,nhbin3,0.,2.,0.) $ ,nhbin2,-0.5,0.5,0.) write(title,504)i,ia call HBOOK1(id_hi_eta4+i*100+ia $ ,title c $ ,nhbin4,0.,3.,0.) c $ ,nhbin3,-1.,1.,0.) $ ,nhbin3,-0.5,.5,0.) enddo enddo return end C............................................................ subroutine fill_eta_histos include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' include '../common/level1.f' include '../common/level2.f' character*10 rzdir * flag to tag low-multilpicity cluster * (according to Samu definition) logical LMSX,LMSY,DRAY 10 format('LADDER',i1) * start loop on planes to fill histos do ip=1,nplanes !1-->6 *--------------------------------------------------------- * main x-y cluster parameters *--------------------------------------------------------- xxx = xm_nt(ip,1) yyy = ym_nt(ip,1) call whichsensor(ip,xxx,yyy,iladder,isensor) if(iladder.eq.0.or.isensor.eq.0)then goto 10102 endif multx = mult(which_cl_x(ip)) multy = mult(which_cl_y(ip)) thetax = axv_nt(ip,1) thetay = ayv_nt(ip,1) LMSX=.false. LMSY=.false. * =============================================== * multiplicity selection of SAMUELE c$$$ if( c$$$ $ (thetax.le.15..and.multx.le.2).or. c$$$ $ (thetax.gt.15..and.multx.le.3).or. c$$$ $ .false.)LMSX=.true. c$$$ if( c$$$ $ (multx.le.2).or. c$$$ $ .false.)LMSY=.true. * =============================================== * looser cut ... after LANDI warning * (I also require that both clx anc cly are low-multiplicity) if( $ (abs(thetax).le.10..and.multx.le.3).or. $ ((abs(thetax).gt.10..and.abs(thetax).le.15.).and.multx.le.4).or. $ (abs(thetax).gt.15..and.multx.le.5).or. $ .false.)LMSX=.true. if( $ (multy.le.3).or. $ .false.)LMSY=.true. * charge scaling factor (to vertical trajectory) pi=acos(-1.) fact= $ sqrt( $ tan(pi*thetax/180.)**2 + $ tan(pi*thetay/180.)**2 + $ 1. $ ) dedxscalx = dedx(which_cl_x(ip)) / fact dedxscaly = dedx(which_cl_y(ip)) / fact dedxscal = (dedxscalx+dedxscaly)/2 DRAY=.false. if(dedxscal.gt.dedxmax)DRAY=.true. *--------------------------------------------------------- write(rzdir,10)iladder call HCDIR('//'//rzdir,'') call HCDIR('//PAWC/'//rzdir,'') *--------------------------------------------------------- *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * X VIEWS *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ *--------------------------------------------------------- if( thetax.lt.angL(1).or. $ thetax.ge.angR(nangbin).or. $ .not.LMSX.or. $ .not.LMSY.or. $ DRAY.or. $ .false.)goto 10101 ! skip iview = (nplanes - ip)*2+2 !12-->2 e2=cog(2,which_cl_x(ip)) e3=cog(3,which_cl_x(ip)) e4=cog(4,which_cl_x(ip)) do ia=1,nangbin if (thetax.lt.angR(ia))goto 111 enddo 111 ibin_ang = ia !angular bin avangbin(ibin_ang,iview)=avangbin(ibin_ang,iview)+thetax navangbin(ibin_ang,iview)=navangbin(ibin_ang,iview)+1 *!!!!!!!!NB non sempre MAXS coincide con la strip massima !!!! if(e2.ge.0.5)e2=e2-1. if(e2.lt.-0.5)e2=e2+1. if(e3.ge.0.5)e3=e3-1. if(e3.lt.-0.5)e3=e3+1. c if(e4.ge.1.)e4=e4-2. c if(e4.lt.-1.)e4=e4+2. if(e4.ge.0.5)e4=e4-1. if(e4.lt.-0.5)e4=e4+1. *!!!!!!!!NB if(e2.ne.0.) $ call hfill(id_hi_eta2+iview*100+ibin_ang, $ e2,0.,1.) if(e3.ne.0.) $ call hfill(id_hi_eta3+iview*100+ibin_ang, $ e3,0.,1.) if(e4.ne.0.) $ call hfill(id_hi_eta4+iview*100+ibin_ang, $ e4,0.,1.) 10101 continue *--------------------------------------------------------- *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * Y VIEWS *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ *--------------------------------------------------------- if( thetay.lt.angL(1).or. $ thetay.ge.angR(nangbin).or. $ .not.LMSY.or. $ .not.LMSX.or. $ DRAY.or. $ .false.)goto 10102 !skip iview = (nplanes - ip)*2+1 !11-->1 e2=cog(2,which_cl_y(ip)) e3=cog(3,which_cl_y(ip)) e4=cog(4,which_cl_y(ip)) do ia=1,nangbin c if (thetay.lt.(ang(ia)+angstep/2.))goto 112 c if (thetay.lt.angbin(ia))goto 112 if (thetay.lt.angR(ia))goto 112 enddo 112 ibin_ang = ia avangbin(ibin_ang,iview)=avangbin(ibin_ang,iview)+thetay navangbin(ibin_ang,iview)=navangbin(ibin_ang,iview)+1 c print*,'****',avangbin(ibin_ang,iview),thetay *!!!!!!!!NB non sempre MAXS coincide con la strip massima !!!! if(e2.gt.0.5)e2=e2-1. if(e2.lt.-0.5)e2=e2+1. if(e3.ge.0.5)e3=e3-1. if(e3.lt.-0.5)e3=e3+1. c if(e4.ge.1.)e4=e4-2. c if(e4.lt.-1.)e4=e4+2. if(e4.ge.0.5)e4=e4-1. if(e4.lt.-0.5)e4=e4+1. *!!!!!!!!NB if(e2.ne.0.) $ call hfill(id_hi_eta2+iview*100+ibin_ang, $ e2,0.,1.) if(e3.ne.0.) $ call hfill(id_hi_eta3+iview*100+ibin_ang, $ e3,0.,1.) if(e4.ne.0.) $ call hfill(id_hi_eta4+iview*100+ibin_ang, $ e4,0.,1.) 10102 continue enddo call HCDIR('//PAWC','') return end C............................................................ subroutine book_charge_histos include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' character*64 title !histos titleinfo -f g77 M LEX c 500 format('Charge correlation (ladder ',i1,' plane ',i2,')') 500 format('Charge correlation ( plane ',i2,')') 501 format('Charge residuals ( plane',i2,')') do ip=1,nplanes write(title,500)ip call HBOOK2(id_hi_chargeco+ip*100+4 $ ,title $ ,200,0.,1000. $ ,200,0.,1000.,0.) write(title,501)ip call HBOOK1(id_hi_chargeco+ip*100 $ ,title $ ,200,-100.,100.,0.) enddo return end C............................................................ subroutine fill_charge_histos(iladder) include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' include '../common/level1.f' include '../common/level2.f' * start loop on planes c parameter (ncorrmax=20000) c real chargexy(nplanes,ncorrmax,2) c integer whichplane,goodtr c common/chargecorr/chargexy,whichplane,goodtr character*10 rzdir 1010 format('LADDER',i1) write(rzdir,1010)iladder call HCDIR('//'//rzdir,'') call HCDIR('//PAWC/'//rzdir,'') do ip=1,nplanes !1-->6 do ie=1,min(goodtr,ncorrmax) sx = chargexy(ip,iladder,ie,1) sy = chargexy(ip,iladder,ie,2) if(sx.ne.0..and.sy.ne.0)then call HFILL(id_hi_chargeco+ip*100+4,sx,sy,1.) d= sy - syksx(ip)*sx - sycsx(ip) d=d/sqrt(syksx(ip)**2+1) call HFILL(id_hi_chargeco+ip*100, $ d,0.,1.) endif enddo enddo return end C............................................................ subroutine book_check_histos * * some histos to check the sample tracks * and define quality cuts * include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' character*64 title !histos title real xbins(nangbin+1) parameter (nrigbin=500) real xrigbin(nrigbin+1) * * * * * * * * * * * * * * * * * * * * * * * * * CHI**2 + RIGIDITY rmin=0.2 rmax=1000. step=abs(log10(rmin)-log10(rmax))/nrigbin xrigbin(1)=rmin do i=2,nrigbin+1 xrigbin(i)=10**(log10(xrigbin(i-1))+step) * print*,'### ',i,xrigbin(i-1),xrigbin(i) enddo call HBOOKB(id_hi_rig $ ,'Rigidity' $ ,nrigbin,xrigbin,0.) call HBOOK1(id_hi_chi $ ,'Chi**2' $ ,1000,0.,100.,0.) call HBPROF(id_hi_chi+100 $ ,'Chi**2 vs Rigidity' $ ,150,rmin,0.1*rmax,0.,100.,0.) call HBPROF(id_hi_chi+200 $ ,'Chi**2 vs Angle x' $ ,50,0.,25.,0.,100.,0.) * * * * * * * * * * * * * * * * * * * * * * * * * ANGLE nbin=200. from=angL(1)!angbin(1) to=angR(nangbin)!angbin(nang+1) do ip=1,nplanes call HBOOK1(id_hi_angle+100+ip $ ,'Angle (Y view)' $ ,nbin,from,to,0.) call HBOOK1(id_hi_angle+200+ip $ ,'Angle (X view)' $ ,nbin,from,to,0.) c$$$ call HBOOK1(id_hi_angle_stat+100+ip c$$$ $ ,'Angle (Y view)' c$$$ $ ,nangbin,from,to,0.) c$$$ call HBOOK1(id_hi_angle_stat+200+ip c$$$ $ ,'Angle (X view)' c$$$ $ ,nangbin,from,to,0.) xbins(1)=angL(1) do i=1,nangbin xbins(i+1)=angR(i) enddo call HBOOKB(id_hi_angle_stat+100+ip $ ,'Angle (Y view)' $ ,nangbin,XBINS,0.) call HBOOKB(id_hi_angle_stat+200+ip $ ,'Angle (X view)' $ ,nangbin,XBINS,0.) enddo * * * * * * * * * * * * * * * * * * * * * * * * * MULTIPLICITY 511 format('Multiplicity (view ',i2,', angle bin ',i2,' )') 512 format('Multiplicity vs angle (view ',i2,' )') do i=1,nviews do ia=1,nangbin write(title,511)i,ia call HBOOK1(id_hi_mult+i*100+ia $ ,title $ ,10,0.5,10.5,0.) enddo write(title,512)i call HBOOK2(id_hi_mult+i*100 $ ,title $ ,50,0.,25.,10,0.5,10.5,0.) call HBOOK2(80000+id_hi_mult+i*100 $ ,title $ ,50,0.,25.,10,0.5,10.5,0.) enddo * * * * * * * * * * * * * * * * * * * * * * * * * AVERAGE CLUSTER 510 format('Average cluster (view ',i2,', angle bin ',i2,' )') do i=1,nviews do ia=1,nangbin write(title,510)i,ia call HBOOK1(id_hi_mult+i*100+ia+50000 $ ,title $ ,int(2*nhside+1) $ ,-(float(nhside)+0.5),(float(nhside)+0.5),0.) enddo enddo * * * * * * * * * * * * * * * * * * * * * * * * * CHARGE 5111 format('dE/dX scaled to 300 micron (view ',i2,' )') 5112 format('dE/dX vs R (view ',i2,' )') 5113 format('dE/dX vs multiplicity (view ',i2,' )') do i=1,nviews write(title,5111)i call HBOOK1(id_hi_charge+i*100 $ ,title $ ,500,0.,1000.,0.) write(title,5112)i call HBPROF(id_hi_charge+i*100+1 $ ,title $ ,1000,0.,100.,0.,1000.,0.) write(title,5113)i call HBOOK2(id_hi_charge+i*100+2 $ ,title $ ,100,0.,1000.,10,0.5,10.5,0.) call HBOOK2(80000+id_hi_charge+i*100+2 $ ,title $ ,100,0.,1000.,10,0.5,10.5,0.) enddo call HBOOK1(id_hi_charge $ ,'dE/dX scaled to 300 micron AVERAGE' $ ,500,0.,1000.,0.) * * * * * * * * * * * * * * * * * * * * * * * * * EVENT SPATIAL DISTRIBUTION 51111 format('Strip Y (plane ',i2,', sensor ',i1,')') 51112 format('Strip X (plane ',i2,', sensor ',i1,')') do ip=1,nplanes do is=1,2 write(title,51111)ip,is call HBOOK1(id_hi_beam+ip*100+is+100000 $ ,title $ ,256,0.5,3072.5,0.) write(title,51112)ip,is call HBOOK1(id_hi_beam+ip*100+is+200000 $ ,title $ ,256,0.5,3072.5,0.) enddo enddo * * * * * * * * * * * * * * * * * * * * * * * * * MAGNETIC FIELD 61 format('Bx (plane ',i2,')') 62 format('By (plane ',i2,')') 63 format('Bz (plane ',i2,')') do ip=1,nplanes write(title,61)ip call HBOOK1(id_hi_bx+ip $ ,title $ ,1000,-5.,5.,0.) write(title,62)ip call HBOOK1(id_hi_by+ip $ ,title $ ,1000,-5.,5.,0.) write(title,63)ip call HBOOK1(id_hi_bz+ip $ ,title $ ,1000,-5.,5.,0.) enddo return end C............................................................ subroutine fill_check_histos include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' include '../common/common_mech.f' include '../common/level1.f' include '../common/level2.f' real v(3),b(3) character*10 rzdir real clusterx(2*nhside+1) real clustery(2*nhside+1) logical DRAY,LMSX,LMSY *--------------------------------------------------------- * track info parameters *--------------------------------------------------------- rig = abs(1./al_nt(5,1)) chi = chi2_nt(1) 1010 format('LADDER',i1) * i plot della rigidita` del chi2 * non dipendono dal ladder... do il=1,3 write(rzdir,1010)il call HCDIR('//'//rzdir,'') call HCDIR('//PAWC/'//rzdir,'') call HFILL(id_hi_chi,chi,0.,1.) call HFILL(id_hi_rig,rig,0.,1.) thetax = axv_nt(1,1) thetay = ayv_nt(1,1) call HFILL(id_hi_chi+100,chi,rig,1.) call HFILL(id_hi_chi+200,chi,thetax,1.) enddo dedxscalav = 0 npdedxscalav = 0 * start loop on planes do ip=1,nplanes !TOP-->BOTTOM *--------------------------------------------------------- * main x-y cluster parameters *--------------------------------------------------------- iviewx = nviewx(7-ip) iviewy = nviewy(7-ip) xxx = xm_nt(ip,1) yyy = ym_nt(ip,1) call whichsensor((7-ip),xxx,yyy,iladder,isensor) if(iladder.eq.0.or.isensor.eq.0)then print*,'*** out *** ',ip,xxx,yyy,iladder,isensor goto 10101 endif multx = mult(which_cl_x(ip)) multy = mult(which_cl_y(ip)) thetax = axv_nt(ip,1) thetay = ayv_nt(ip,1) * =============================================== * looser multiplicity cut ... after LANDI warning * (I also require that both clx anc cly are low-multiplicity) LMSX=.false. LMSY=.false. if( $ (abs(thetax).le.10..and.multx.le.3).or. $ ((abs(thetax).gt.10..and.abs(thetax).le.15.).and.multx.le.4).or. $ (abs(thetax).gt.15..and.multx.le.5).or. $ .false.)LMSX=.true. if( $ (multy.le.3).or. $ .false.)LMSY=.true. * charge scaling factor (to vertical trajectory) pi=acos(-1.) fact= $ sqrt( $ tan(pi*thetax/180.)**2 + $ tan(pi*thetay/180.)**2 + $ 1. $ ) dedxscalx = dedx(which_cl_x(ip)) / fact dedxscaly = dedx(which_cl_y(ip)) / fact dedxscal = (dedxscalx+dedxscaly)/2 DRAY=.false. if(dedxscal.gt.dedxmax)DRAY=.true. if(.not.DRAY.and.LMSX.and.LMSX)then dedxscalav = dedxscalav + dedxscal npdedxscalav = npdedxscalav + 1 endif ista = indstart(which_cl_x(ip)) ice = indmax(which_cl_x(ip)) ilast = totCLlength if(which_cl_x(ip).ne.nclstr1) $ ilast = indstart(which_cl_x(ip)+1)-1 do i = -nhside,nhside ii = ice + i iii=i+nhside+1 clusterx(iii) = 0 if(ii.le.ilast.and.ii.ge.ista) $ clusterx(iii) = clsignal(ii) / dedx(which_cl_x(ip)) enddo ista = indstart(which_cl_y(ip)) ice = indmax(which_cl_y(ip)) ilast = totCLlength if(which_cl_y(ip).ne.nclstr1) $ ilast = indstart(which_cl_y(ip)+1)-1 do i = -nhside,nhside ii = ice + i iii=i+nhside+1 clustery(iii) = 0 if(ii.le.ilast.and.ii.ge.ista) $ clustery(iii) = clsignal(ii) / dedx(which_cl_y(ip)) enddo *--------------------------------------------------------- write(rzdir,1010)iladder call HCDIR('//'//rzdir,'') call HCDIR('//PAWC/'//rzdir,'') *--------------------------------------------------------- * find the angular bins (ibin_angy ibin_angx) * associated to the track intersection point with plane ip *--------------------------------------------------------- ia=0 if( thetax.lt.angL(1).or. $ thetax.ge.angR(nangbin))goto 111 do ia=1,nangbin if (thetay.lt.angR(ia))goto 111 enddo 111 ibin_angy = ia ia=0 if( thetay.lt.angL(1).or. $ thetay.ge.angR(nangbin))goto 112 do ia=1,nangbin if (thetax.lt.angR(ia))goto 112 enddo 112 ibin_angx = ia *--------------------------------------------------------- * fill histos *--------------------------------------------------------- ************************************************ ANGLE call HFILL(id_hi_angle+200+ip,thetax,0.,1.) call HFILL(id_hi_angle+100+ip,thetay,0.,1.) call HFILL(id_hi_angle_stat+200+ip,thetax,0.,1.) call HFILL(id_hi_angle_stat+100+ip,thetay,0.,1.) ************************************************ multiplicity & cluster if(ibin_angx.ne.0)then call HFILL(id_hi_mult+iviewx*100+ibin_angx, $ float(mult(which_cl_x(ip))),0.,1.) do i=-nhside,nhside ii = i + nhside +1 xxxx=clusterx(ii) call HFILL(id_hi_mult+iviewx*100+ibin_angx+50000, $ float(i),0.,xxxx) enddo endif if(ibin_angy.ne.0)then call HFILL(id_hi_mult+iviewy*100+ibin_angy, $ float(mult(which_cl_y(ip))),0.,1.) do i=-nhside,nhside ii = i + nhside +1 xxxx=clustery(ii) call HFILL(id_hi_mult+iviewy*100+ibin_angy+50000, $ float(i),0.,xxxx) enddo endif ************************************************ charge call HFILL(id_hi_charge+iviewx*100, $ dedxscalx,0.,1.) call HFILL(id_hi_charge+iviewy*100, $ dedxscaly,0.,1.) call HFILL(id_hi_charge+iviewx*100+1, $ rig,dedxscalx,1.) call HFILL(id_hi_charge+iviewy*100+1, $ rig,dedxscaly,1.) * -------------------------------------------- C this are the condition applied to select the C particle sample for pfa calibration * -------------------------------------------- if(.not.DRAY.and.LMSX.and.LMSY)then call HFILL(id_hi_charge+iviewx*100+2, $ dedxscalx,float(multx),1.) call HFILL(id_hi_charge+iviewy*100+2, $ dedxscaly,float(multy),1.) call HFILL(id_hi_mult+iviewx*100, $ abs(thetax),float(multx),1.) call HFILL(id_hi_mult+iviewy*100, $ abs(thetay),float(multy),1.) endif call HFILL(80000+id_hi_charge+iviewx*100+2, $ dedxscalx,float(multx),1.) call HFILL(80000+id_hi_charge+iviewy*100+2, $ dedxscaly,float(multy),1.) call HFILL(80000+id_hi_mult+iviewx*100, $ abs(thetax),float(multx),1.) call HFILL(80000+id_hi_mult+iviewy*100, $ abs(thetay),float(multy),1.) ************************************************ beam profle call hfill(id_hi_beam+ip*100+isensor+100000, $ float(MAXS(which_cl_y(ip))),0.,1.) call hfill(id_hi_beam+ip*100+isensor+200000, $ float(MAXS(which_cl_x(ip))),0.,1.) ************************************************ magnetic field v(1) = xv_nt(ip,1) v(2) = yv_nt(ip,1) v(3) = zv_nt(ip,1) call gufld(v,b) bx=b(1) by=b(2) bz=b(3) call hfill(id_hi_bx+ip,bx,0.,1.) call hfill(id_hi_by+ip,by,0.,1.) call hfill(id_hi_bz+ip,bz,0.,1.) enddo dedxscalav = dedxscalav/npdedxscalav call HFILL(id_hi_charge, $ dedxscalav,0.,1.) 10101 continue call HCDIR('//PAWC','') return end *************************************************************************** * * * * * * * * * *************************************************************************** subroutine FCNcharge(npar,grad,fval,xval,iflag,futil) include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' double precision grad(*),xval(*),fval integer npar,iflag double precision KKK,CCC c parameter (ncorrmax=20000) c real chargexy(nplanes,ncorrmax,2) c integer whichplane,goodtr c common/chargecorr/chargexy,whichplane,goodtr KKK=xval(1) CCC=xval(2) fval = 0. sum2=0 do i=1,min(goodtr,ncorrmax) sx=chargexy(whichplane,whichladder,i,1) sy=chargexy(whichplane,whichladder,i,2) * if(sx.gt.140..and.sy.gt.140)then d=(sy-KKK*sx-CCC)**2 d=d/(KKK**2+1) d=d/(1.9*4.**2+1.2*8.**2) !errors sum2=sum2+d c$$$ PRINT*,'xx ',i,d,' sx-sy ' c$$$ $ ,chargexy(whichplane,whichladder,i,1) c$$$ $ ,chargexy(whichplane,whichladder,i,2) * endif enddo fval = sum2 return end ************************************************** ************************************************** * CHARGE CORRELATION CALIBRATION ************************************************** ************************************************** subroutine charge_calibration(iladder) include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' c include '../common/common_mech.f' c include '../common/level1.f' c include '../common/level2.f' external FCNcharge character*12 chnam double precision val,error do ip=1,nplanes print*,' ' print*,' ' print*,' ' print*, $ ' -----------------------------------' print*, $ '////////////// MINUIT for plane ',ip,' - ladder ',iladder print*, $ ' -----------------------------------' open(UNIT=40,FILE= $ 'bin-aux/data-card.dat' $ ,STATUS='OLD') call mintio(40,6,7) whichplane = ip whichladder = iladder call minuit(FCNcharge,0) c call MNERRS(1,eplus,eminus,eparab,globcc) call MNPOUT(1,chnam,val,error,bnd1,bnd2,ivarbl) syksx(ip)=val esyksx(ip)=error call MNPOUT(2,chnam,val,error,bnd1,bnd2,ivarbl) sycsx(ip)=val esycsx(ip)=error close(40) enddo return end ************************************************** ************************************************** * ETA CALIBRATION ************************************************** ************************************************** subroutine eta_calibration(iladder,ia) include '../common/commontracker.f' include '../common/calib.f' include '../common/common_calibration.f' * -------------------------------------------- * some local variables for eta correction real histo2(nhbin2) c real eta2corr(nhbin2+1,nviews) !eta2 correction c real eta2abs(nhbin2+1) !eta2 abscissa * -------------------------------------------- character*64 ctemp character*10 rzdir 1010 format('LADDER',i1) write(rzdir,1010)iladder call HCDIR('//'//rzdir,'') call HCDIR('//PAWC/'//rzdir,'') do i=1,nviews id=id_hi_eta2+i*100+ia * eta2 histogram call HUNPAK(id,histo2,' ',0) !eta2 histogram call HGIVE(id,ctemp,nx,xmi,xma,ny,ymi,yma,nwt,loc) sum=HSUM(id) c print*,ctemp,nx,xmi,xma,ny,ymi,yma,nwt,loc if(i.eq.1)then do ie=0,nx eta2abs(ie+1)=xmi+ie*(xma-xmi)/nx enddo endif if(sum.eq.0)then do ie=0,nx eta2corr(ie+1,i)=eta2abs(ie+1) enddo goto 3 endif * ------------------------------------------------------ * if the histogram contains less than NMIN4ETA2 entries, * the correction factor is not calculated * CHANGED!!! * it is calculated anyway.... * ------------------------------------------------------ nmin4eta2 = 200 if(sum.lt.nmin4eta2)then c$$$ do ie=0,nx c$$$ eta2corr(ie+1,i)=eta2abs(ie+1) c$$$ enddo print*,'*** WARNING ***' print*,'ETA2 parameters: less than ',nmin4eta2 $ ,' entries ----> ',rzdir,' view ',i,' ang. bin ',ia c$$$ goto 3 endif * integral of the eta2 histogram eta2corr(1,i)=0. do ie=1,nx eta2corr(ie+1,i)=eta2corr(ie,i)+histo2(ie) enddo * eta2 normalization angle=avangbin(ia,i) x0=eta2fix(i,angle) !point of reference for correction function do ifix=1,nx+1 if(eta2abs(ifix).ge.x0)goto 2 enddo 2 y0=eta2corr(ifix,i)/sum c shift=eta2shift(i,angle) shift = 0. do ie=0,nx eta2corr(ie+1,i)= $ x0-y0+eta2corr(ie+1,i)/sum-shift enddo 3 continue enddo !end loop on views return end *---***---***---***---***---***---***---***---*** * * functions to define the point where to "fix" * the eta correction * and the shift * *---***---***---***---***---***---***---***---*** function eta2fix(iview,angle) if(mod(iview,2).eq.0)then !>>>>>>>> X eta2fix = -0.5 else !>>>>>>>> Y eta2fix = -0.5 endif return end c$$$ function eta2shift(iview,angle) c$$$ c$$$ parameter (hallm_e=1650.) !cm**2 / Vs c$$$ parameter (hallm_h=310.) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_preanalysis.f' c$$$ character*32 ctemp c$$$ real hallm c$$$ c$$$ pi=acos(-1.) c$$$ c$$$ np=npl(iview) c$$$ if(mod(iview,2).eq.0)then !>>>>>>>> X c$$$ id=id_hi_by+np c$$$ hallm=hallm_h*1e-4 c$$$ else c$$$ id=id_hi_bx+np c$$$ hallm=hallm_e*1e-4 c$$$ endif c$$$ bfield=HSTATI(id,1,' ',0) c$$$ bfield=bfield*1e-1 c$$$ tgth_h=-hallm*bfield c$$$ c$$$ tgth_eff=tan(pi*angle/180.)-tgth_h c$$$ theta_eff=atan(tgth_eff)*180./pi c$$$ if(iview.ne.12)theta_eff=-theta_eff c$$$ c$$$ if(mod(iview,2).eq.0)then !>>>>>>>> X c$$$ c$$$ eta2shift = -1.549*sin(2*pi*theta_eff/13.59) c$$$ eta2shift = eta2shift/51. c$$$ c$$$ else !>>>>>>>> Y c$$$ c$$$ eta2shift = -1.348*sin(2*pi*theta_eff/27.73) c$$$ eta2shift = eta2shift/67. c$$$* NBNBNBNBNB ... da risolvere il problema della simulazione!!!! c$$$ c$$$ endif c$$$ c$$$c print*,iview,np,bfield,angle,theta_eff,'>>>>>>> ',eta2shift c$$$ c$$$ return c$$$ end