Bug Summary

File:libraries/TAGGER/DTAGMHit_factory_Calib.cc
Location:line 190, column 13
Description:Value stored to 'pulse_peak' is never read

Annotated Source Code

1// $Id$
2//
3// File: DTAGMHit_factory_Calib.cc
4// Created: Sat Aug 2 12:23:43 EDT 2014
5// Creator: jonesrt (on Linux gluex.phys.uconn.edu)
6//
7
8
9#include <iostream>
10#include <iomanip>
11using namespace std;
12
13#include "DTAGMDigiHit.h"
14#include "DTAGMTDCDigiHit.h"
15#include "DTAGMGeometry.h"
16#include "DTAGMHit_factory_Calib.h"
17#include <DAQ/Df250PulseIntegral.h>
18#include <DAQ/Df250PulsePedestal.h>
19#include <DAQ/Df250Config.h>
20
21using namespace jana;
22
23const int DTAGMHit_factory_Calib::k_fiber_good;
24const int DTAGMHit_factory_Calib::k_fiber_bad;
25const int DTAGMHit_factory_Calib::k_fiber_noisy;
26
27//------------------
28// init
29//------------------
30jerror_t DTAGMHit_factory_Calib::init(void)
31{
32 DELTA_T_ADC_TDC_MAX = 10.0; // ns
33 USE_ADC = 0;
34 CUT_FACTOR = 1;
35 gPARMS->SetDefaultParameter("TAGMHit:DELTA_T_ADC_TDC_MAX", DELTA_T_ADC_TDC_MAX,
36 "Maximum difference in ns between a (calibrated) fADC time and"
37 " F1TDC time for them to be matched in a single hit");
38 gPARMS->SetDefaultParameter("TAGMHit:CUT_FACTOR", CUT_FACTOR, "TAGM pulse integral cut factor, 0 = no cut");
39 gPARMS->SetDefaultParameter("TAGMHit:USE_ADC", USE_ADC, "Use ADC times in TAGM");
40
41 CHECK_FADC_ERRORS = false;
42 gPARMS->SetDefaultParameter("TAGMHit:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits");
43
44 // initialize calibration constants
45 fadc_a_scale = 0;
46 fadc_t_scale = 0;
47 t_base = 0;
48 t_tdc_base=0;
49
50 // calibration constants stored in row, column format
51 for (int row = 0; row <= TAGM_MAX_ROW5; ++row) {
52 for (int col = 0; col <= TAGM_MAX_COLUMN102; ++col) {
53 fadc_gains[row][col] = 0;
54 fadc_pedestals[row][col] = 0;
55 fadc_time_offsets[row][col] = 0;
56 tdc_time_offsets[row][col] = 0;
57 fiber_quality[row][col] = 0;
58 }
59 }
60
61 return NOERROR;
62}
63
64//------------------
65// brun
66//------------------
67jerror_t DTAGMHit_factory_Calib::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
68{
69 // Only print messages for one thread whenever run number change
70 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
71 static set<int> runs_announced;
72 pthread_mutex_lock(&print_mutex);
73 bool print_messages = false;
74 if(runs_announced.find(runnumber) == runs_announced.end()){
75 print_messages = true;
76 runs_announced.insert(runnumber);
77 }
78 pthread_mutex_unlock(&print_mutex);
79
80 /// set the base conversion scales
81 fadc_a_scale = 1.1; // pixels per count
82 fadc_t_scale = 0.0625; // ns per count
83 t_base = 0.; // ns
84
85 pthread_mutex_unlock(&print_mutex);
86 if(print_messages) jout << "In DTAGMHit_factory_Calib, loading constants..." << std::endl;
87
88 // load base time offset
89 map<string,double> base_time_offset;
90 if (eventLoop->GetCalib("/PHOTON_BEAM/microscope/base_time_offset",base_time_offset))
91 jout << "Error loading /PHOTON_BEAM/microscope/base_time_offset !" << endl;
92 if (base_time_offset.find("TAGM_BASE_TIME_OFFSET") != base_time_offset.end())
93 t_base = base_time_offset["TAGM_BASE_TIME_OFFSET"];
94 else
95 jerr << "Unable to get TAGM_BASE_TIME_OFFSET from /PHOTON_BEAM/microscope/base_time_offset !" << endl;
96 if (base_time_offset.find("TAGM_TDC_BASE_TIME_OFFSET") != base_time_offset.end())
97 t_tdc_base = base_time_offset["TAGM_TDC_BASE_TIME_OFFSET"];
98 else
99 jerr << "Unable to get TAGM_TDC_BASE_TIME_OFFSET from /PHOTON_BEAM/microscope/base_time_offset !" << endl;
100
101 if (load_ccdb_constants("fadc_gains", "gain", fadc_gains) &&
102 load_ccdb_constants("fadc_pedestals", "pedestal", fadc_pedestals) &&
103 load_ccdb_constants("fadc_time_offsets", "offset", fadc_time_offsets) &&
104 load_ccdb_constants("tdc_time_offsets", "offset", tdc_time_offsets) &&
105 load_ccdb_constants("fiber_quality", "code", fiber_quality) &&
106 load_ccdb_constants("tdc_timewalk_corrections", "c0", tw_c0) &&
107 load_ccdb_constants("tdc_timewalk_corrections", "c1", tw_c1) &&
108 load_ccdb_constants("tdc_timewalk_corrections", "c2", tw_c2) &&
109 load_ccdb_constants("tdc_timewalk_corrections", "threshold", tw_c3) &&
110 load_ccdb_constants("tdc_timewalk_corrections", "reference", ref) &&
111 load_ccdb_constants("integral_cuts", "integral", int_cuts))
112 {
113 return NOERROR;
114 }
115 return UNRECOVERABLE_ERROR;
116}
117
118//------------------
119// evnt
120//------------------
121jerror_t DTAGMHit_factory_Calib::evnt(JEventLoop *loop, uint64_t eventnumber)
122{
123 /// Generate DTAGMHit object for each DTAGMDigiHit object.
124 /// This is where the first set of calibration constants
125 /// is applied to convert from digitzed units into natural
126 /// units.
127 ///
128 /// Note that this code does NOT get called for simulated
129 /// data in HDDM format. The HDDM event source will copy
130 /// the precalibrated values directly into the _data vector.
131
132 // extract the TAGM geometry
133 vector<const DTAGMGeometry*> tagmGeomVect;
134 eventLoop->Get( tagmGeomVect );
135 if (tagmGeomVect.size() < 1)
136 return OBJECT_NOT_AVAILABLE;
137 const DTAGMGeometry& tagmGeom = *(tagmGeomVect[0]);
138
139 const DTTabUtilities* locTTabUtilities = nullptr;
140 loop->GetSingle(locTTabUtilities);
141
142 // First loop over all TAGMDigiHits and make DTAGMHits out of them
143 vector<const DTAGMDigiHit*> digihits;
144 loop->Get(digihits);
145 for (unsigned int i=0; i < digihits.size(); i++) {
146 const DTAGMDigiHit *digihit = digihits[i];
147
148 // Throw away hits with firmware errors (post-summer 2016 firmware)
149 if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF))
150 continue;
151
152 // Get pedestal, prefer associated event pedestal if it exists,
153 // otherwise, use the average pedestal from CCDB
154 double pedestal = fadc_pedestals[digihit->row][digihit->column];
155 double nsamples_integral = (double)digihit->nsamples_integral;
156 double nsamples_pedestal = (double)digihit->nsamples_pedestal;
157
158 // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
159 if(nsamples_pedestal == 0) {
160 jerr << "DTAGMDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl;
161 continue;
162 }
163
164 if ( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) {
165 // the measured pedestal must be scaled by the ratio of the number
166 // of samples used to calculate the integral and the pedestal
167 // Changed to conform to D. Lawrence changes Dec. 4 2014
168 double ped_sum = (double)digihit->pedestal;
169 pedestal = ped_sum * nsamples_integral/nsamples_pedestal;
170 }
171 double single_sample_ped = pedestal/nsamples_pedestal;
172
173 double pulse_peak = 0.0;
174 if(digihit->datasource == 1) { // handle pre-Fall 2016 firmware
175 // Throw away hits where the fADC timing algorithm failed
176 //if (digihit->pulse_time == 0) continue;
177 // The following condition signals an error state in the flash algorithm
178 // Do not make hits out of these
179 const Df250PulsePedestal* PPobj = nullptr;
180 digihit->GetSingle(PPobj);
181 if (PPobj != nullptr){
182 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
183 pulse_peak = PPobj->pulse_peak - PPobj->pedestal;
184 }
185
186 // Skip events where fADC algorithm fails
187 //if (digihit->pulse_time == 0) continue; // Should already be caught above
188 } else {
189 // starting with the Fall 2016 firmware, we can get all of the values directly from the digihit
190 pulse_peak = digihit->pulse_peak - single_sample_ped;
Value stored to 'pulse_peak' is never read
191 }
192
193 // throw away hits from bad or noisy fibers
194 int quality = fiber_quality[digihit->row][digihit->column];
195 if (quality == k_fiber_bad || quality == k_fiber_noisy)
196 continue;
197
198 // Skip digihit if pulse peak is lower than cut value
199 double P = digihit->pulse_peak - pedestal/nsamples_pedestal;
200
201 // Apply calibration constants
202 double A = digihit->pulse_integral;
203 double T = digihit->pulse_time;
204 A -= pedestal;
205 int row = digihit->row;
206 int column = digihit->column;
207 if (A < CUT_FACTOR*int_cuts[row][column]) continue;
208
209 DTAGMHit *hit = new DTAGMHit;
210 hit->row = row;
211 hit->column = column;
212 double Elow = tagmGeom.getElow(column);
213 double Ehigh = tagmGeom.getEhigh(column);
214 hit->E = (Elow + Ehigh)/2;
215
216 hit->integral = A;
217 hit->pulse_peak = P;
218 hit->npix_fadc = A * fadc_a_scale * fadc_gains[row][column];
219 hit->time_fadc = T * fadc_t_scale - fadc_time_offsets[row][column] + t_base;
220 hit->t=hit->time_fadc;
221 hit->has_TDC=false;
222 hit->has_fADC=true;
223
224 hit->AddAssociatedObject(digihit);
225 _data.push_back(hit);
226 }
227
228 // Next, loop over TDC hits, matching them to the existing fADC hits
229 // where possible and updating their time information. If no match is
230 // found, then create a new hit with just the TDC info.
231 vector<const DTAGMTDCDigiHit*> tdcdigihits;
232 loop->Get(tdcdigihits);
233 for (unsigned int i=0; i < tdcdigihits.size(); i++) {
234 const DTAGMTDCDigiHit *digihit = tdcdigihits[i];
235
236 // Apply calibration constants here
237 int row = digihit->row;
238 int column = digihit->column;
239 double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - tdc_time_offsets[row][column] + t_tdc_base;
240
241 // Look for existing hits to see if there is a match
242 // or create new one if there is no match
243 DTAGMHit *hit = nullptr;
244 for (unsigned int j=0; j < _data.size(); ++j) {
245 if (_data[j]->row == row && _data[j]->column == column &&
246 fabs(T - _data[j]->time_fadc) < DELTA_T_ADC_TDC_MAX)
247 {
248 hit = _data[j];
249 }
250 }
251 if (hit == nullptr) {
252 hit = new DTAGMHit;
253 hit->row = row;
254 hit->column = column;
255 double Elow = tagmGeom.getElow(column);
256 double Ehigh = tagmGeom.getEhigh(column);
257 hit->E = (Elow + Ehigh)/2;
258 hit->time_fadc = 0;
259 hit->npix_fadc = 0;
260 hit->integral = 0;
261 hit->pulse_peak = 0;
262 hit->has_fADC=false;
263 _data.push_back(hit);
264 }
265 hit->time_tdc=T;
266 hit->has_TDC=true;
267
268 // apply time-walk corrections
269 double P = hit->pulse_peak;
270 double c0 = tw_c0[row][column];
271 double c1 = tw_c1[row][column];
272 double c2 = tw_c2[row][column];
273 //double TH = thresh[row][column];
274 double c3 = tw_c3[row][column];
275 double t0 = ref[row][column];
276 //pp_0 = TH*pow((pp_0-c0)/c1,1/c2);
277 if (P > 0) {
278 //T -= c1*(pow(P/TH,c2)-pow(pp_0/TH,c2));
279 T -= c1*pow(1/(P+c3),c2) - (t0 - c0);
280 }
281
282 // Optionally only use ADC times
283 if (USE_ADC && hit->has_fADC) hit->t = hit->time_fadc;
284 else hit->t = T;
285
286 hit->AddAssociatedObject(digihit);
287 }
288
289 return NOERROR;
290}
291
292//------------------
293// erun
294//------------------
295jerror_t DTAGMHit_factory_Calib::erun(void)
296{
297 return NOERROR;
298}
299
300//------------------
301// fini
302//------------------
303jerror_t DTAGMHit_factory_Calib::fini(void)
304{
305 return NOERROR;
306}
307
308//---------------------
309// load_ccdb_constants
310//---------------------
311bool DTAGMHit_factory_Calib::load_ccdb_constants(
312 std::string table_name,
313 std::string column_name,
314 double result[TAGM_MAX_ROW5+1][TAGM_MAX_COLUMN102+1])
315{
316 std::vector< std::map<std::string, double> > table;
317 std::string ccdb_key = "/PHOTON_BEAM/microscope/" + table_name;
318 if (eventLoop->GetCalib(ccdb_key, table))
319 {
320 jout << "Error loading " << ccdb_key << " from ccdb!" << std::endl;
321 return false;
322 }
323 for (unsigned int i=0; i < table.size(); ++i) {
324 int row = (table[i])["row"];
325 int col = (table[i])["column"];
326 result[row][col] = (table[i])[column_name];
327 }
328 return true;
329}