Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JRateK40.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <map>
4
5#include "TROOT.h"
6#include "TFile.h"
7#include "TH1D.h"
8#include "TFormula.h"
9
10#include "JROOT/JRootToolkit.hh"
11
13
14#include "JPhysics/KM3NeT.hh"
15#include "JPhysics/KM3NeT2D.hh"
16#include "JPhysics/Antares.hh"
17#include "JPhysics/JPDF.hh" //for module radius
19
20#include "Jeep/JContainer.hh"
21#include "Jeep/JParser.hh"
22#include "Jeep/JMessage.hh"
23
24
25/**
26 * \file
27 *
28 * Example program to calculate singles rate.\n
29 * The calculation is based on Antares internal note ANTARES-PHYS-2006-005 by Juergen Brunner.\n
30 * According Antares internal note ANTARES-PHYS-2012-013, the absorption length is used (and not the attenuation length).
31 * \author mdejong, vkulikovskiy
32 */
33int main(int argc, char* argv[])
34{
35 using namespace std;
36 using namespace JPP;
37
38 typedef JContainer< map<double, double> > container_type;
39
40 string outputFile;
41 double bequerel;
42 JPMTParameters parameters;
43 container_type fd;
44 TFormula f1;
45 int debug;
46
47 try {
48
49 JProperties properties = parameters.getProperties();
50
51 JParser<> zap("Example program to calculate singles rate.");
52
53 zap['o'] = make_field(outputFile) = "k40.root";
54 zap['b'] = make_field(bequerel, "radioactivity") = 13.75e3; // [m^-3 s^-1]
55 zap['P'] = make_field(properties, "PMT parameters") = JPARSER::initialised();
56 zap['f'] = make_field(fd, "filter data") = JPARSER::initialised();
57 zap['F'] = make_field(f1, "filter function") = JPARSER::initialised();
58 zap['d'] = make_field(debug) = 3;
59
60 zap(argc, argv);
61 }
62 catch(const exception &error) {
63 FATAL(error.what() << endl);
64 }
65
66
67 using namespace NAMESPACE;
68
70
71 for (const auto& i : fd) {
72 g1[i.first] = i.second;
73 }
74
75 g1.compile();
76 g1.setExceptionHandler(new typename JPolint1Function1D_t::JDefaultResult(JMATH::zero));
77
78 const JPMTAnalogueSignalProcessor cpu(parameters);
79
80 //For Geant4 inputs see https://git.km3net.de/vkulikovskiy/g4ly
81
82 const double wmin = 280.0; // minimal wavelength [nm]
83 const double wmax = 700.0; // maximal wavelength [nm]
84 const double ng = 47.11; // Geant4 simulations of K40 decays in (wmin,wmax) window
85 const int npe = 1; // number of photo-electrons for each decay
86 const double cpow = 2.156; // 1/lambda^cpow shape of the Cherenkov photons distribution is a fit from Geant4
87
88 TFile out(outputFile.c_str(), "recreate");
89
90 const double dx = 1.5; // [nm]
91 const int nx = (int) ((wmax - wmin) / dx);
92
93 TH1D h0("h0", NULL, nx, wmin, wmax);
94 TH1D h1("h1", NULL, nx, wmin, wmax);
95
96
97 double R[] = { 0.0, 0.0 };
98 double W[] = { 0.0, 0.0 };
99
100 const string option[] = { "1Dx1D", "2D" };
101
102 for (double x = -1.0, dx = 0.02; x <= +1.0; x += dx) {
103 W[0] += getPhotocathodeArea() * getAngularAcceptance(x) * dx;
104 }
105
106 double Y = 1.0;
107
108 Y *= cpow - 1.0;
109 Y /= pow(wmin, 1.0 - cpow) - pow(wmax, 1.0 - cpow);
110
111 Y *= bequerel * ng;
112
113 Y *= cpu.getSurvivalProbability(npe);
114 Y *= 0.5e-3;
115
116
117 double A[] = { 0.0, 0.0 };
118
119 for (int ix = 1; ix <= h0.GetXaxis()->GetNbins(); ++ix) {
120
121 const double w = h0.GetXaxis()->GetBinCenter(ix);
122 const double dw = h0.GetXaxis()->GetBinWidth (ix);
123
124 W[1] = 0.0;
125
126 for (double x = -1.0, dx = 0.02; x <= +1.0; x += dx) {
127 W[1] += KM3NET2D::getPhotocathodeArea2D(x, w) * dx;
128 }
129
130 double U = Y / pow(w,cpow);
131
132 if (!g1.empty()) {
133 U *= g1(w);
134 }
135 if (f1.IsValid()) {
136 U *= f1.Eval(w);
137 }
138
139 R[0] += W[0] * U * dw * getQE(w) * getAbsorptionLength(w);
140 R[1] += W[1] * U * dw * getAbsorptionLength(w);
141 A[0] += W[0] * U * dw * getQE(w);
142 A[1] += W[1] * U * dw;
143
144 h0.SetBinContent(ix, W[0] * U * getQE(w) * getAbsorptionLength(w));
145 h1.SetBinContent(ix, W[1] * U * getAbsorptionLength(w));
146 }
147
148 double labs[] = { 0.0, 0.0 };
149
150 cout << "PMT survival probability: " << cpu.getSurvivalProbability(npe) << endl;
151
152 for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
153
154 labs[i] = R[i]/A[i];
155
156 const double corr = R[i]*(1.0-exp(-MODULE_RADIUS_M/labs[i])); //rate due to coincidences within the DOM
157
158 cout << setw(6) << left << option[i] + ":" << right
159 << " rate " << FIXED(7,3) << R[i] << " [kHz]"
160 << " corrected " << FIXED(7,3) << R[i] - corr << " [kHz]"
161 << " <l_abs> " << FIXED(7,3) << labs[i] << " [m]" << endl;
162 }
163
164 out.Write();
165 out.Close();
166}
Properties of Antares PMT and deep-sea water.
Container I/O.
string outputFile
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition JDrawLED.cc:68
double getAbsorptionLength(const double lambda)
Definition JDrawPD0.cc:27
General purpose messaging.
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
PMT analogue signal processor.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Double_t g1(const Double_t x)
Function.
Definition JQuantiles.cc:25
int main(int argc, char *argv[])
Definition JRateK40.cc:33
Properties of KM3NeT PMT and deep-sea water.
Properties of KM3NeT PMT and deep-sea water.
Data structure for PMT parameters.
JProperties getProperties(const JEquationParameters &equation=JPMTParameters::getEquationParameters())
Get properties of this class.
Utility class to parse parameter values.
Utility class to parse command line options.
Definition JParser.hh:1698
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for effective photo-cathode area of PMT.
Definition KM3NeT2D.hh:5235
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
virtual double getSurvivalProbability(const int NPE) const override
Probability that a hit survives the simulation of the PMT.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
void setExceptionHandler(const JSupervisor &supervisor)
Set the supervisor for handling of exceptions.
Type definition of a 1st degree polynomial interpolation with result type double.