New default nplotter_ktopanom_new.f90

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