Merge branch 'main' into issue_156

Author: c-dilks

bin/environ.sh

@@ -42,7 +42,7 @@ groovy_classpath=(
 timeline_java_opts=(
   -DCLAS12DIR=$COATJAVA/
   -Djava.util.logging.config.file=$COATJAVA/etc/logging/$log_config.properties
-  -Xmx1024m
+  -Xmx1536m
   -XX:+UseSerialGC
 )
 timeline_groovy_opts=(

detectors/pom.xml

@@ -51,7 +51,7 @@
     <dependency>
       <groupId>org.apache.groovy</groupId>
       <artifactId>groovy-all</artifactId>
-      <version>4.0.25</version>
+      <version>4.0.26</version>
       <type>pom</type>
       <exclusions>
         <exclusion>
@@ -64,7 +64,7 @@
     <dependency>
       <groupId>org.apache.groovy</groupId>
       <artifactId>groovy-dateutil</artifactId>
-      <version>4.0.25</version>
+      <version>4.0.26</version>
     </dependency>
     <!-- https://mvnrepository.com/artifact/org.codehaus.gpars/gpars -->
     <dependency>

detectors/src/main/java/org/jlab/clas/timeline/util/HistoUtil.groovy

@@ -65,25 +65,46 @@ class HistoUtil {
     if( !(hist.getEntries() > 0) || !(hist.integral() > 0) ) { // `.getEntries()` can be nonzero for empty histograms, so be sure to check the integral too
       return defIfEmpty
     }
-    def hist_list = []
-    hist.getAxis().getNBins().times { bin ->
-      def counts = hist.getBinContent(bin).toInteger() // FIXME: assumes the histogram is unweighted
+    def bins_list = []
+    def cum_counts_list = []
+    int sum = 0
+    int nbins = hist.getAxis().getNBins()
+    nbins.times { bin ->
+      int counts = hist.getBinContent(bin).toInteger() // FIXME: assumes the histogram is unweighted
       def value  = hist.getAxis().getBinCenter(bin)
-      counts.times { hist_list += value }
+      bins_list += value
+      sum += counts
+      cum_counts_list += sum
     }
-    def listMedian = { d ->
-      if(d.size()==0) {
+    def listMedian = { _sum, bins, cum_counts ->
+      if(bins.size()==0 || cum_counts.size()==0) {
         // this list may end up being empty if there are _few_ entries in the histogram, since this method
         // is called to get the quartiles; in this case, just return `defIfEmpty`
         return defIfEmpty
       }
-      d.sort()
-      def m = d.size().intdiv(2)
-      d.size() % 2 ? d[m] : (d[m-1]+d[m]) / 2
+      // Compute the median using the fact that you have histogrammed data to do this slightly more efficiently
+      int mid_count = _sum.intdiv(2)
+      boolean is_even = (_sum%2==0)
+      if (is_even) mid_count += 1 //NOTE: If you have an even length list there is no middle element and you will need to average the middle+1 and middle-1 elements.
+      for (int i=0; i<cum_counts.size(); i++) {
+        if (cum_counts[i]+1==mid_count && is_even) {
+          return [(bins[i]+bins[i+1])/2, i+1, mid_count]
+        }
+        if (cum_counts[i]>mid_count) { 
+          return [bins[i], i, mid_count]
+        }
+      }
+      return [bins[bins.size()-1], bins.size()-1, mid_count] //NOTE: SHOULD NEVER REACH THIS POINT
     }
-    def mq = listMedian(hist_list)
-    def lq = listMedian(hist_list.findAll{it<mq})
-    def uq = listMedian(hist_list.findAll{it>mq})
+
+    def mq_list = listMedian(sum,bins_list,cum_counts_list)
+    int mq_idx = mq_list[1]
+    int mq_sum = mq_list[2]
+    def lq_list = listMedian(mq_sum,bins_list.subList(0,mq_idx+1),cum_counts_list.subList(0,mq_idx+1))
+    def lq = lq_list[0]
+    def uq_list = listMedian(mq_sum*3,bins_list.subList(mq_idx,nbins),cum_counts_list.subList(mq_idx,nbins))
+    def uq = uq_list[0]
+
     return uq - lq
   }
 }

qa-physics/monitorPlot.groovy

@@ -143,15 +143,19 @@ def buildMonAveGr = { tObj ->
   return gr
 }
 
-// build monitor 'average X' distribution
-def buildMonAveDist = { tObj,nb,lb,ub ->
+// build monitor 'runSum' distribution (sum of bins' distribtions for a run)
+def buildMonRunSumDist = { tObj ->
   def histN = objToMonName(tObj.getName())
   def histT = objToMonTitle(tObj.getTitle())
-  histN = histN + "_aveDist"
-  histT = "average " + histT
+  histN = histN + "_runSumDist"
   histT = histT.replaceAll(/$/,' distribution')
   histT = appendLegend(histT)
-  def hist = new H1F(histN,histT,nb,lb,ub)
+  def hist = new H1F(
+    histN,
+    histT,
+    tObj.getAxis().getNBins(),
+    tObj.getAxis().min(),
+    tObj.getAxis().max())
   if(histN.contains("_hm_")) { hist.setLineColor(2); }
   return hist
 }
@@ -206,8 +210,6 @@ def objList
 def part
 def hel
 def varStr
-def varNB
-def varLB,varUB
 def runnum
 def timeBinNum
 def obj
@@ -245,8 +247,8 @@ inList.each { inFile ->
       T.addLeaf(monTree,[runnum,'helic','sinPhi',part,hel,'aveGr'],{
         buildMonAveGr(obj)
       })
-      T.addLeaf(monTree,[runnum,'helic','sinPhi',part,hel,'aveDist'],{
-        buildMonAveDist(obj,100,-0.25,0.25)
+      T.addLeaf(monTree,[runnum,'helic','sinPhi',part,hel,'runSumDist'],{
+        buildMonRunSumDist(obj)
       })
 
       // instantiate sinPhi dists binned for an asymmetry, denoted "asymGrid" (if not
@@ -267,7 +269,7 @@ inList.each { inFile ->
         aveXerr = obj.getRMS() / Math.sqrt(ent)
         monTree[runnum]['helic']['sinPhi'][part][hel]['aveGr'].addPoint(
           timeBinNum, aveX, 0, aveXerr )
-        monTree[runnum]['helic']['sinPhi'][part][hel]['aveDist'].fill(aveX)
+        monTree[runnum]['helic']['sinPhi'][part][hel]['runSumDist'].add(obj)
         // add rebinned <sinPhi> distribution to asymGrid
         obj.getAxis().getNBins().times { bin ->
           def counts = obj.getBinContent(bin)
@@ -292,10 +294,10 @@ inList.each { inFile ->
         return g
       })
       T.addLeaf(monTree,[runnum,'helic','dist','heldef','heldefDist'],{
-        def h = buildMonAveDist(obj,50,0,1)
-        def hN = h.getName().replaceAll(/_aveDist$/,'_heldefDist')
-        h.setName(hN)
-        h.setTitle('average defined helicity fraction distribution')
+        def histN = objToMonName(obj.getName()) + "_heldefDist"
+        def histT = appendLegend('average defined helicity fraction distribution')
+        def h = new H1F(histN, histT, 50, 0, 1)
+        if(histN.contains("_hm_")) { h.setLineColor(2); }
         return h
       })
       ent = obj.integral()
@@ -324,10 +326,10 @@ inList.each { inFile ->
         return g
       })
       T.addLeaf(monTree,[runnum,'helic','dist','rellum','rellumDist'],{
-        def h = buildMonAveDist(obj,50,0.9,1.1)
-        def hN = h.getName().replaceAll(/_aveDist$/,'_rellumDist')
-        h.setName(hN)
-        h.setTitle('average n+/n- distribution')
+        def histN = objToMonName(obj.getName()) + "_rellumDist"
+        def histT = appendLegend('average n+/n- distribution')
+        def h = new H1F(histN, histT, 50, 0.9, 1.1)
+        if(histN.contains("_hm_")) { h.setLineColor(2); }
         return h
       })
       if(obj.integral()>0) {
@@ -379,13 +381,7 @@ inList.each { inFile ->
         varStr = tok[1]
 
         T.addLeaf(monTree,[runnum,'DIS',varStr,'aveGr'],{buildMonAveGr(obj)})
-        T.addLeaf(monTree,[runnum,'DIS',varStr,'aveDist'],{
-          varNB = 100
-          if(varStr=="Q2") { varLB=0; varUB=12; }
-          else if(varStr=="W") { varLB=0; varUB=6; }
-          else { varLB=0; varUB=1; }
-          buildMonAveDist(obj,varNB,varLB,varUB)
-        })
+        T.addLeaf(monTree,[runnum,'DIS',varStr,'runSumDist'],{buildMonRunSumDist(obj)})
 
         // add <varStr> point to the monitors
         ent = obj.integral()
@@ -394,7 +390,7 @@ inList.each { inFile ->
           aveXerr = obj.getRMS() / Math.sqrt(ent)
           monTree[runnum]['DIS'][varStr]['aveGr'].addPoint(
             timeBinNum, aveX, 0, aveXerr )
-          monTree[runnum]['DIS'][varStr]['aveDist'].fill(aveX)
+          monTree[runnum]['DIS'][varStr]['runSumDist'].add(obj)
         }
       }
     }
@@ -405,24 +401,14 @@ inList.each { inFile ->
       part = tok[1]
       varStr = tok[2]
       T.addLeaf(monTree,[runnum,'inclusive',part,varStr,'aveGr'],{buildMonAveGr(obj)})
-      T.addLeaf(monTree,[runnum,'inclusive',part,varStr,'aveDist'],{
-        varNB = 100
-        varLB = 0
-        varUB = 0
-        if(varStr=='p') { varLB=0; varUB=10 }
-        else if(varStr=='pT') { varLB=0; varUB=4 }
-        else if(varStr=='z') { varLB=0; varUB=1 }
-        else if(varStr=='theta') { varLB=0; varUB=Math.toRadians(90.0) }
-        else if(varStr=='phiH') { varLB=-3.15; varUB=3.15 }
-        buildMonAveDist(obj,varNB,varLB,varUB)
-      })
+      T.addLeaf(monTree,[runnum,'inclusive',part,varStr,'runSumDist'],{buildMonRunSumDist(obj)})
       ent = obj.integral()
       if(ent>0) {
         aveX = obj.getMean()
         aveXerr = obj.getRMS() / Math.sqrt(ent)
         monTree[runnum]['inclusive'][part][varStr]['aveGr'].addPoint(
           timeBinNum, aveX, 0, aveXerr )
-        monTree[runnum]['inclusive'][part][varStr]['aveDist'].fill(aveX)
+        monTree[runnum]['inclusive'][part][varStr]['runSumDist'].add(obj)
       }
     }
 
@@ -470,7 +456,7 @@ inList.each { inFile ->
 // build timelines
 //---------------------
 
-// loop through 'aveDist' monitors: for each one, add its mean to the timeline
+// loop through monitors: for each one, add its mean to the timeline
 def timelineTree = [:]
 T.exeLeaves(monTree,{
   if(T.key.contains('Dist') || T.key.contains('Graph')) {
@@ -552,7 +538,7 @@ T.exeLeaves(monTree,{
     }
 
     // add this run's <X> to the timeline (and stddev to the stddev timelines)
-    if(T.key=='aveDist') {
+    if(T.key=='runSumDist') {
       aveX = T.leaf.getMean()
       stddevX = T.leaf.getRMS()
       aveXerr = stddevX / Math.sqrt(T.leaf.integral())