DetectorSensitivity/Functors.h at master · SkyToGround/DetectorSensitivity · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
//
//  Functors.h
//

#ifndef Functors_h
#define Functors_h

#include "Detector.h"

struct FindActivityFunctor : DenseFunctor<double> {
	ArrayXd S;
	double B;
	unsigned int C_L;
	double beta;
	FindActivityFunctor(ArrayXd S, double B, unsigned int C_L, double beta) : DenseFunctor<double> (1, 1), S(S), B(B), C_L(C_L), beta(beta) {

	}

	//Should yield the same result as CalcTruePositiveProb
	int operator()(const VectorXd &p, VectorXd &res) const {
		double M = p[0];
		ArrayXd tot = S * M + B;
		res = VectorXd(1);
		ArrayXd tgt = ArrayXd::Zero(tot.size());
		if (C_L > 100) {
			ArrayXi i = ArrayXi::LinSpaced(C_L - 1, 1, C_L - 1);
			for (int j = 0; j < tot.size(); j++) {
				tgt[j] = exp(-tot[j] + i.cast<double>()*log(tot[j]) - log((1.0 + 1.0/(12.0*i.cast<double>()) + 1.0/(288.0*i.cast<double>()*i.cast<double>()))*sqrt(2*pi*i.cast<double>()))- i.cast<double>()*log(i.cast<double>()/e)).sum();
			}
		} else {
			ArrayXi i = ArrayXi::LinSpaced(C_L, 0, C_L - 1);
			for (int j = 0; j < tot.size(); j++) {
				tgt[j] = exp(-tot[j])*(pow(tot[j], i.cast<double>())/factorial(i)).sum();
			}
		}
		res[0] = tgt.prod() - beta;
		return 0;
	}
};

struct FindActivityFunctorLM : DenseFunctor<double> {
	Detector *det;
	unsigned int C_L;
	double beta;
	unsigned int iterations;
	FindActivityFunctorLM(Detector *det, unsigned int C_L, double beta, unsigned int iterations) : DenseFunctor<double> (1, 1), det(det), C_L(C_L), beta(beta), iterations(iterations) {

	}

	int operator()(const VectorXd &p, VectorXd &res) const {
		double M = p[0];
		res[0] = det->SimMeasurementsOpti(M, C_L, iterations) - beta;
		return 0;
	}
};

#endif /* Functors_h */