ctlike.cpp
4.04 KB
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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
/***************************************************************************
* ctlike - CTA maximum likelihood tool *
* ----------------------------------------------------------------------- *
* copyright (C) 2010 by Jurgen Knodlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
/**
* @file ctlike.cpp
* @brief CTA maximum likelihood tool implementation
* @author J. Knodlseder
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdio.h>
#include "ctlike.hpp"
#include "GTools.hpp"
/***********************************************************************//**
* @brief Run maximum likelihood analysis
***************************************************************************/
int ctlike::run(void)
{
// Initialise return code
int rc = 0;
// Get standard parameters
m_method = toupper(par("method")->value());
m_caldb = par("caldb")->value();
m_irf = par("irf")->value();
m_srcmdl = par("srcmdl")->value();
m_models = GModels(m_srcmdl);
// Set fixed parameters
m_max_iter = 1000;
m_opt = new GOptimizerLM;
// Set optimizer parameters
((GOptimizerLM*)m_opt)->max_iter(m_max_iter);
// Branch on method
if (m_method == "UNBINNED")
rc = unbinned();
else
rc = binned();
// Free optimizer
delete m_opt;
// Return
return rc;
}
/***********************************************************************//**
* @brief Perform unbinned maximum likelihood analysis
***************************************************************************/
int ctlike::unbinned(void)
{
// Initialise return code
int rc = 0;
// Get parameters
std::string evfile = par("evfile")->value();
// Declare observations
GObservations obs;
GCTAObservation run;
// Setup ROI covered by data
GCTAInstDir instDir;
GCTARoi roi;
instDir.radec_deg(117.02, -33.35); // Adapt to file
roi.centre(instDir);
roi.radius(2.5);
// Setup energy range covered by data
GEnergy emin;
GEnergy emax;
emin.TeV(0.02);
emax.TeV(100.0);
// Setup time range covered by data
GTime tstart;
GTime tstop;
tstart.met(0.0);
tstop.met(10000.0);
// Load data and response and set ROI, energy range and time range
// for analysis
run.load_unbinned(evfile);
run.response(m_irf, m_caldb);
run.roi(roi);
run.gti()->add(tstart, tstop);
run.ebounds()->append(emin, emax);
obs.append(run);
// Setup models for optimizing
obs.models(m_models);
// Perform LM optimization
obs.optimize(*m_opt);
std::cout << obs << std::endl;
// Return
return rc;
}
/***********************************************************************//**
* @brief Perform binned maximum likelihood analysis
***************************************************************************/
int ctlike::binned(void)
{
// Initialise return code
int rc = 0;
// Get parameters
std::string cntmap = par("cntmap")->value();
// Declare observations
GObservations obs;
GCTAObservation run;
// Load binned CTA observation
run.load_binned(cntmap);
run.response(m_irf, m_caldb);
obs.append(run);
// Setup models for optimizing
obs.models(m_models);
// Perform LM optimization
obs.optimize(*m_opt);
std::cout << obs << std::endl;
// Return
return rc;
}