1 package org.lcsim.analysis.util;
2
3 import hep.aida.ICloud1D;
4 import hep.aida.IEvaluator;
5 import hep.aida.ITuple;
6 import java.util.Arrays;
7 import org.lcsim.util.aida.AIDA;
8
9
10
11
12
13 public class BetterRMS90Calculator {
14
15 public Result calculateRMS90(ICloud1D cloud) {
16 return calculate(cloud2array(cloud));
17 }
18
19 public Result calculateRMS90(ITuple tuple, String expression) {
20 AIDA aida = AIDA.defaultInstance();
21 IEvaluator eval = aida.analysisFactory().createTupleFactory(aida.tree()).createEvaluator(expression);
22 return calculate(tuple2array(tuple, eval));
23 }
24
25 public Result calculateRMS90(ITuple tuple, IEvaluator evaluator) {
26 return calculate(tuple2array(tuple, evaluator));
27 }
28
29 private double[] tuple2array(ITuple tuple, IEvaluator evaluator) {
30 evaluator.initialize(tuple);
31 int entries = tuple.rows();
32 double[] elist = new double[entries];
33 tuple.start();
34 for (int i = 0; i < entries; i++) {
35 tuple.next();
36 elist[i] = evaluator.evaluateDouble();
37 }
38 return elist;
39 }
40
41 private double[] cloud2array(ICloud1D cloud) {
42 int entries = cloud.entries();
43 double[] elist = new double[entries];
44 for (int i = 0; i < entries; i++) {
45 elist[i] = cloud.value(i);
46 }
47 return elist;
48 }
49
50 Result calculate(double[] elist) {
51 int entries = elist.length;
52 double rms90 = Double.NaN;
53 double mean90 = Double.NaN;
54 int ntail = (int) (.1 * entries);
55 int ncore = entries - ntail;
56 Arrays.sort(elist);
57
58
59 double svi = 0.;
60 double sv2i = 0.;
61 for (int i = ntail; i < ncore; i++) {
62 svi += elist[i];
63 sv2i += elist[i] * elist[i];
64 }
65
66 for (int k = 0; k <= ntail; k++) {
67 double sv = svi;
68 double sv2 = sv2i;
69 for (int i = k; i < ntail; i++) {
70 sv += elist[i];
71 sv2 += elist[i] * elist[i];
72 }
73 for (int i = ncore; i < ncore+k; i++) {
74 sv += elist[i];
75 sv2 += elist[i] * elist[i];
76 }
77 double mean = sv / ncore;
78 double rms = Math.sqrt(sv2 / ncore - mean * mean);
79 if (Double.isNaN(rms90) || rms < rms90) {
80 rms90 = rms;
81 mean90 = mean;
82 }
83 }
84 return new Result(rms90, mean90);
85 }
86
87
88 }