1 | /* |
---|
2 | * TimePixCalibrationHandler.cpp |
---|
3 | * |
---|
4 | * Created on: March, 2012 |
---|
5 | * Author: John Idarraga <idarraga@cern.ch> |
---|
6 | */ |
---|
7 | |
---|
8 | |
---|
9 | #include "TimePixCalibrationHandler.h" |
---|
10 | |
---|
11 | #include <iostream> |
---|
12 | #include <fstream> |
---|
13 | |
---|
14 | #include <stdlib.h> |
---|
15 | |
---|
16 | using namespace std; |
---|
17 | |
---|
18 | namespace MAFTools { |
---|
19 | |
---|
20 | TimePixCalibrationHandler::TimePixCalibrationHandler(const char * af, |
---|
21 | const char * bf, |
---|
22 | const char * cf, |
---|
23 | const char * tf, |
---|
24 | int width, int height){ |
---|
25 | |
---|
26 | m_width = width; |
---|
27 | m_height = height; |
---|
28 | |
---|
29 | ifstream afs(af, fstream::in); |
---|
30 | ifstream bfs(bf, fstream::in); |
---|
31 | ifstream cfs(cf, fstream::in); |
---|
32 | ifstream tfs(tf, fstream::in); |
---|
33 | |
---|
34 | //////////////////////////////////// |
---|
35 | // a |
---|
36 | m_a = new double[m_width*m_height]; |
---|
37 | double temp = 0.; |
---|
38 | int i = 0; |
---|
39 | while ( afs.good() ) { |
---|
40 | afs >> temp; |
---|
41 | if(afs.eof()) break; |
---|
42 | m_a[i++] = temp; |
---|
43 | } |
---|
44 | afs.close(); |
---|
45 | if(i != m_width*m_height){ |
---|
46 | cout << "Calibration does not seem to be complete for paremeter a." << endl; |
---|
47 | cout << "Got " << i << " items. " << m_width << "*" << m_height |
---|
48 | << " were requested. Giving up." << endl; |
---|
49 | exit(1); |
---|
50 | } |
---|
51 | |
---|
52 | //////////////////////////////////// |
---|
53 | // b |
---|
54 | m_b = new double[m_width*m_height]; |
---|
55 | i = 0; |
---|
56 | while ( bfs.good() ) { |
---|
57 | bfs >> temp; |
---|
58 | if(bfs.eof()) break; |
---|
59 | m_b[i++] = temp; |
---|
60 | } |
---|
61 | bfs.close(); |
---|
62 | if(i != m_width*m_height){ |
---|
63 | cout << "Calibration does not seem to be complete for paremeter b." << endl; |
---|
64 | cout << "Got " << i << " items. " << m_width << "*" << m_height |
---|
65 | << " were requested. Giving up." << endl; |
---|
66 | exit(1); |
---|
67 | } |
---|
68 | |
---|
69 | //////////////////////////////////// |
---|
70 | // c |
---|
71 | m_c = new double[m_width*m_height]; |
---|
72 | i = 0; |
---|
73 | while ( cfs.good() ) { |
---|
74 | cfs >> temp; |
---|
75 | if(cfs.eof()) break; |
---|
76 | m_c[i++] = temp; |
---|
77 | } |
---|
78 | cfs.close(); |
---|
79 | if(i != m_width*m_height){ |
---|
80 | cout << "Calibration does not seem to be complete for paremeter c." << endl; |
---|
81 | cout << "Got " << i << " items. " << m_width << "*" << m_height |
---|
82 | << " were requested. Giving up." << endl; |
---|
83 | exit(1); |
---|
84 | } |
---|
85 | |
---|
86 | //////////////////////////////////// |
---|
87 | // t |
---|
88 | m_t = new double[m_width*m_height]; |
---|
89 | i = 0; |
---|
90 | while ( tfs.good() ) { |
---|
91 | tfs >> temp; |
---|
92 | if(tfs.eof()) break; |
---|
93 | m_t[i++] = temp; |
---|
94 | } |
---|
95 | tfs.close(); |
---|
96 | if(i != m_width*m_height){ |
---|
97 | cout << "Calibration does not seem to be complete for paremeter t." << endl; |
---|
98 | cout << "Got " << i << " items. " << m_width << "*" << m_height |
---|
99 | << " were requested. Giving up." << endl; |
---|
100 | exit(1); |
---|
101 | } |
---|
102 | |
---|
103 | } |
---|
104 | |
---|
105 | TimePixCalibrationHandler::~TimePixCalibrationHandler(){ |
---|
106 | } |
---|
107 | |
---|
108 | double TimePixCalibrationHandler::GetE(pair<int,int> pix, int tot){ |
---|
109 | |
---|
110 | //TF1 * surr = GetSurrogateTF1(pix); |
---|
111 | |
---|
112 | int index = XYtoC(pix, m_width); |
---|
113 | |
---|
114 | double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] }; |
---|
115 | |
---|
116 | /* |
---|
117 | if(pix.first == 111 && pix.second == 70) { |
---|
118 | std::cout << "a = " << par[0] << std::endl; |
---|
119 | std::cout << "b = " << par[1] << std::endl; |
---|
120 | std::cout << "c = " << par[2] << std::endl; |
---|
121 | std::cout << "t = " << par[3] << std::endl; |
---|
122 | }*/ |
---|
123 | |
---|
124 | |
---|
125 | double evalTOT = 0.; |
---|
126 | double prev_e = 0.; |
---|
127 | double dist = 0.; |
---|
128 | double prev_dist = 10000; |
---|
129 | |
---|
130 | // search |
---|
131 | for(double e = 3.0; e <= 10000. ; e+=0.05) { // steps of 0.05 keV, start search |
---|
132 | |
---|
133 | //evalTOT = surr->Eval(e); |
---|
134 | evalTOT = SurrogateCalibFunction(e, par); |
---|
135 | if(evalTOT < 0) continue; // don't consider those values. Not physical. |
---|
136 | dist = TMath::Abs(evalTOT - tot); |
---|
137 | //if(pix.first == 111 && pix.second == 70) cout << "dist = " << dist << " | e = " << e << " | evalTOT = " << evalTOT << " | tot = " << tot << endl; |
---|
138 | if( dist > prev_dist ) break; // if I am going far, break. |
---|
139 | |
---|
140 | prev_dist = dist; |
---|
141 | prev_e = e; // keep the last energy |
---|
142 | } |
---|
143 | |
---|
144 | return prev_e; |
---|
145 | } |
---|
146 | |
---|
147 | int TimePixCalibrationHandler::GetTOT (pair<int, int> pix, double E) { |
---|
148 | |
---|
149 | int index = XYtoC(pix, m_width); |
---|
150 | double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] }; |
---|
151 | |
---|
152 | return SurrogateCalibFunction(E, par); |
---|
153 | //return GetSurrogateTF1(pix)->Eval(E); |
---|
154 | } |
---|
155 | |
---|
156 | double TimePixCalibrationHandler::SurrogateCalibFunction(double x, double * par){ |
---|
157 | return par[0]*x + par[1] - (par[2]/(x-par[3])); |
---|
158 | } |
---|
159 | |
---|
160 | #ifndef __MAFALDA_EMBEDDED |
---|
161 | // Only if needed one can fetch a few functions as TF1 objects |
---|
162 | // Not included for embedded applications |
---|
163 | TF1 * TimePixCalibrationHandler::GetSurrogateTF1(pair<int, int> pix){ |
---|
164 | |
---|
165 | // check first if the function has been defined. Use the Set. |
---|
166 | m_surrogateItr = m_surrogateSet.find(pix); |
---|
167 | |
---|
168 | // If surrogate TF1 does not exist yet, build the function first |
---|
169 | if( m_surrogateItr == m_surrogateSet.end() ) { |
---|
170 | TString ftitle = "surrogate_"; |
---|
171 | ftitle += pix.first; ftitle += "_"; ftitle += pix.second; |
---|
172 | // Usually the parameters come given in keV !!! WARNING !!! |
---|
173 | m_surrogateMap[pix] = new TF1(ftitle,"[0]*x + [1] - ([2]/(x-[3]))",0,1000); // 1 MeV |
---|
174 | // Populate parameters. Look for right index first |
---|
175 | int index = XYtoC(pix, m_width); |
---|
176 | m_surrogateMap[pix]->SetParameters(m_a[index], m_b[index], m_c[index], m_t[index]); |
---|
177 | } |
---|
178 | |
---|
179 | return m_surrogateMap[pix]; |
---|
180 | } |
---|
181 | #endif |
---|
182 | |
---|
183 | // These two functions are already defined in MAFTools but for certain tests |
---|
184 | // I want this class to be self-contained. |
---|
185 | int TimePixCalibrationHandler::XYtoC(pair<int, int> position, int dimX) { |
---|
186 | return XYtoC(position.first, position.second, dimX); |
---|
187 | } |
---|
188 | int TimePixCalibrationHandler::XYtoC(int x, int y, int dimX){ |
---|
189 | return y * dimX + x; |
---|
190 | } |
---|
191 | |
---|
192 | } // end namespace MAFTools |
---|