diff --git a/src/User/nplotter_ktopanom_new.f90 b/src/User/nplotter_ktopanom_new.f90 index 463e87f..c3d0281 100644 --- a/src/User/nplotter_ktopanom_new.f90 +++ b/src/User/nplotter_ktopanom_new.f90 @@ -22,9 +22,12 @@ module nplotter_ktopanom implicit none include 'mpicommon.f' - allocate(histos(1)) + allocate(histos(4)) histos(1) = plot_setup_uniform(0.0_dp, 200.0_dp, 5.0_dp, 'mbl') + histos(2) = plot_setup_uniform(0._dp,300._dp,10._dp,'ptl') + histos(3) = plot_setup_uniform(0._dp,300._dp,5._dp,'ptb1') + histos(4) = plot_setup_uniform(0._dp,1._dp,1._dp,'totfid') end subroutine @@ -55,26 +58,35 @@ module nplotter_ktopanom integer :: bjetindices(3) integer :: nonbjetindices(3) - integer :: bjetindices_bbfrac(3) - integer :: nonbjetindices_bbfrac(3) - real(dp) :: plep(4), pnu(4) integer :: i,n integer :: imax, imaxb - real(dp) :: puremass - real(dp) :: mbl + real(dp) :: puremass, ptpure + real(dp) :: mbl, ptb1 + + integer :: numb, numnonb jetindices(:) = [5,6,7] + ! initialize with zero + vals = [0._dp,0._dp,0._dp,0._dp] + wts = [0._dp,0._dp,0._dp,0._dp] + ids = histos + ! the 'bbpart' variable (logical, true or false) designates whether we ! have two b-quarks in the final state or not bjetindices = 0 nonbjetindices = 0 + ! require exactly two jets, return with default zero + if (jets /= 2) then + return + endif + ! for main part with at most one b if (bbpart .eqv. .false.) then imax = 1 @@ -104,11 +116,11 @@ module nplotter_ktopanom endif - ! require at least one b-jet and one light jet - if ( bjetindices(1) == 0 .or. nonbjetindices(1) == 0 ) then - ids = histos - vals = [0._dp] - wts = [0._dp] + numb = count(bjetindices(1:jets) /= 0) + numnonb = count(nonbjetindices(1:jets) /= 0) + + ! require exactly one b-jet and one light jet, return with default zero + if ((numb /= 1) .or. (numnonb /= 1)) then return endif @@ -120,11 +132,13 @@ module nplotter_ktopanom pnu(:) = p(3,:) endif + ! take b-index as leading b-jet + ptb1 = ptpure(p(bjetindices(1),:)) mbl = puremass(p(bjetindices(1),:) + plep(:)) ids = histos - vals = [mbl] - wts = [wt] + vals = [mbl, ptpure(plep), ptb1, 0.5_dp] + wts = [wt, wt, wt, wt] end subroutine book