Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DTACHit_factory.cc
Go to the documentation of this file.
1 /*
2  * DTACHIT_factory.cc
3  *
4  * Created on: Mar 24, 2017
5  * Author: hovanes
6  */
7 /*
8  * DTACHIT_factory.cc
9  *
10  * Created on: Mar 24, 2017
11  * Author: Hovanes Egiyan
12  */
13 
14 #include <cmath>
15 using namespace std;
16 
17 #include "DTACHit_factory.h"
18 
19 // The parameter below can be overwritten from the command line
20 // Define the static data members
22 // Flag to use time-walk corrections
24 // TDC time window in ns
26 // Maximum difference between ADc and TDC time to form a single hit
28 
29 // End JANA command line parameter values
30 
31 // List of announced runs, starts empty
33 
34 //------------------
35 // init
36 //------------------
37 jerror_t DTACHit_factory::init(void) {
38  // Set default configuration parameters
39 
40 // gPARMS->SetDefaultParameter("TAC:CHECK_FADC_ERRORS", checkErrorsOnFADC,
41 // "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits");
42 //
43 // gPARMS->SetDefaultParameter("TAC:HIT_TIME_WINDOW", timeWindowTDC,
44 // "Time window of trigger corrected TDC time in which a hit in"
45 // " in the TDC will match to a hit in the fADC to form an TAC hit");
46 //
47 // gPARMS->SetDefaultParameter("TAC:DELTA_T_ADC_TDC_MAX",
48 // timeDifferencInADCandTDC,
49 // "Maximum difference in ns between a (calibrated) fADC time and"
50 // " CAEN TDC time for them to be matched in a single hit");
51 //
52 // gPARMS->SetDefaultParameter("TAC:USE_TIMEWALK_CORRECTION",
53 // useTimeWalkCorrections,
54 // "Flag to decide if time-walk corrections should be applied.");
55 
56  cout << "In DTACHit_factory::init" << endl;
57  return NOERROR;
58 }
59 
60 //------------------
61 // brun
62 //------------------
63 jerror_t DTACHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) {
64 
65  // Only print messages for one thread whenever run number change
66  bool printMessages = false;
67  if (addRun(runnumber)) {
68  printMessages = true;
69  }
70  /// Read in calibration constants
71  if (printMessages)
72  jout << "In DTACHit_factory, loading constants..." << endl;
73 
74  readCCDB( eventLoop );
75 
76  return NOERROR;
77 }
78 
79 //------------------
80 // evnt
81 //------------------
82 jerror_t DTACHit_factory::evnt(jana::JEventLoop *eventLoop, uint64_t eventNumber) {
83 // cout << "Building basic DTACHit objects" << endl;
84  makeFADCHits(eventLoop, eventNumber);
85  makeTDCHits(eventLoop, eventNumber);
86 
87  return NOERROR;
88 }
89 
90 jerror_t DTACHit_factory::readCCDB(jana::JEventLoop *eventLoop) {
91  // load scale factors
92  map<string, double> scaleFactors;
93  // a_scale (TAC_ADC_SCALE)
94  if (eventLoop->GetCalib("/TAC/digi_scales", scaleFactors))
95  jout << "Error loading /TAC/digi_scales !" << endl;
96  if (scaleFactors.find("TAC_ADC_ASCALE") != scaleFactors.end())
97  energyScale = scaleFactors["TAC_ADC_ASCALE"];
98  else
99  jerr << "Unable to get TAC_ADC_ASCALE from /TAC/digi_scales !" << endl;
100  // t_scale (TAC_ADC_SCALE)
101  if (scaleFactors.find("TAC_ADC_TSCALE") != scaleFactors.end())
102  timeScaleADC = scaleFactors["TAC_ADC_TSCALE"];
103  else
104  jerr << "Unable to get TAC_ADC_TSCALE from /TAC/digi_scales !" << endl;
105 
106  // load base time offset
107  map<string, double> baseTimeOffsets;
108  // t_base (TAC_BASE_TIME_OFFSET)
109  if (eventLoop->GetCalib("/TAC/base_time_offset", baseTimeOffsets))
110  jout << "Error loading /TAC/base_time_offset !" << endl;
111  if (baseTimeOffsets.find("TAC_BASE_TIME_OFFSET") != baseTimeOffsets.end())
112  timeBaseADC = baseTimeOffsets["TAC_BASE_TIME_OFFSET"];
113  else
114  jerr
115  << "Unable to get TAC_BASE_TIME_OFFSET from /TAC/base_time_offset !"
116  << endl;
117  // t_tdc_base (TAC_TDC_BASE_TIME_OFFSET)
118  if (baseTimeOffsets.find("TAC_TDC_BASE_TIME_OFFSET")
119  != baseTimeOffsets.end())
120  timeBaseTDC = baseTimeOffsets["TAC_TDC_BASE_TIME_OFFSET"];
121  else
122  jerr
123  << "Unable to get TAC_TDC_BASE_TIME_OFFSET from /TAC/base_time_offset !"
124  << endl;
125 
126  // load constant tables
127  // a_gains (gains)
128  if (eventLoop->GetCalib("/TAC/gains", energyGain))
129  jout << "Error loading /TAC/gains !" << endl;
130  // a_pedestals (pedestals)
131  if (eventLoop->GetCalib("/TAC/pedestals", adcPedestal))
132  jout << "Error loading /TAC/pedestals !" << endl;
133  // adc_time_offsets (adc_timing_offsets)
134  if (eventLoop->GetCalib("/TAC/adc_timing_offsets", adcTimeOffset))
135  jout << "Error loading /TAC/adc_timing_offsets !" << endl;
136  // tdc_time_offsets (tdc_timing_offsets)
137  if (eventLoop->GetCalib("/TAC/timing_offsets", tdcTimeOffsets))
138  jout << "Error loading /TAC/timing_offsets !" << endl;
139  // timewalk_parameters (timewalk_parms)
140  if (eventLoop->GetCalib("TAC/timewalk_parms", timeWalkParameters))
141  jout << "Error loading /TAC/timewalk_parms !" << endl;
142 
143  return NOERROR;
144 }
145 
146 void DTACHit_factory::makeFADCHits(jana::JEventLoop *eventLoop,
147  uint64_t eventNumber) {
148  const DTTabUtilities* locTTabUtilities = nullptr;
149  eventLoop->GetSingle(locTTabUtilities);
150 
151  // Get TAC hits
152  vector<const DTACDigiHit*> tacDigiHits;
153  eventLoop->Get(tacDigiHits);
154 
155  for (auto& tacDigiHit : tacDigiHits) {
156  // Throw away hits with firmware errors (post-summer 2016 firmware)
157  if (errorCheckIsNeededForFADC()
158  && !locTTabUtilities->CheckFADC250_NoErrors(
159  tacDigiHit->getQF()))
160  continue;
161 
162  // Initialize pedestal to one found in CCDB, but override it
163  // with one found in event if is available (?)
164  // For now, only keep events with a correct pedestal
165  double pedestal = adcPedestal;
166  double nsamples_integral = tacDigiHit->getNsamplesIntegral();
167  double nsamples_pedestal = tacDigiHit->getNsamplesPedestal();
168 
169  // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
170  if (nsamples_pedestal == 0) {
171  jerr << "DSCDigiHit with nsamples_pedestal == 0 ! Event = "
172  << eventNumber << endl;
173  continue;
174  }
175 
176  // tacDigiHit->pedestal is the sum of "nsamples_pedestal" samples
177  // Calculate the average pedestal per sample
178  if ((tacDigiHit->getPedestal() > 0)
179  && locTTabUtilities->CheckFADC250_PedestalOK(tacDigiHit->QF)
180  && nsamples_pedestal > 0) {
181  pedestal = (double) tacDigiHit->getPedestal() / nsamples_pedestal;
182  }
183 
184  // Subtract pedestal from pulse peak
185  if (tacDigiHit->getPulseTime() == 0 || tacDigiHit->getPedestal() == 0
186  || tacDigiHit->getPulsePeak() == 0)
187  continue;
188  double pulse_peak = tacDigiHit->getPulsePeak() - pedestal;
189 
190  // Subtract pedestal from pulse integral
191  double pulseIntegral = (double) tacDigiHit->getPulseIntegral()
192  - pedestal * nsamples_integral;
193 
194  double pulseTime = (double) tacDigiHit->getPulseTime();
195 
196  DTACHit *tacHit = new DTACHit;
197 
198  tacHit->setE(energyScale * energyGain * pulseIntegral);
199 
200  tacHit->setTimeFADC(
201  timeScaleADC * pulseTime - adcTimeOffset + timeBaseADC);
202  tacHit->setT(tacHit->getTimeFADC());
203  tacHit->setTimeTDC(numeric_limits<double>::quiet_NaN());
204 
205  tacHit->setTDCPresent(false);
206  tacHit->setFADCPresent(true);
207 
208  tacHit->setTimeTDC(tacHit->getTimeFADC()); // set time from fADC in case no TDC tacHit
209  tacHit->setPulsePeak(pulse_peak);
210 
211 // vector<pair<string,string>> outputStrings;
212 // tacHit->toStrings( outputStrings );
213 // cout << endl << "Newly build TAC hit " << endl;
214 // for( auto& outputPair : outputStrings ) {
215 // cout << outputPair.first << " = " << outputPair.second << endl;
216 // }
217 
218  tacHit->AddAssociatedObject(tacDigiHit);
219 
220  _data.push_back(tacHit);
221  }
222 
223  return;
224 }
225 
226 void DTACHit_factory::makeTDCHits(jana::JEventLoop *eventLoop, uint64_t eventnumber) {
227  // Next, loop over TDC hits, matching them to the
228  // existing fADC hits where possible and updating
229  // their time information. If no match is found, then
230  // create a new hit with just the TDC info.
231  const DTTabUtilities* locTTabUtilities = nullptr;
232  eventLoop->GetSingle(locTTabUtilities);
233 
234  vector<const DTACTDCDigiHit*> tacTDCDigiHits;
235  eventLoop->Get(tacTDCDigiHits);
236 
237  for (auto& tacTDCDigiHit : tacTDCDigiHits) {
238 
239  // Get the the time and corrected it for the TDC base
240  double tdcTime = locTTabUtilities->Convert_DigiTimeToNs_CAEN1290TDC(
241  tacTDCDigiHit) - timeBaseTDC;
242 
243  // Look for existing hits to see if there is a match
244  // or create new one if there is no match
245  // Require that the trigger corrected TDC time fall within
246  // a reasonable time window so that when a hit is associated with
247  // a hit in the TDC and not the ADC it is a "decent" TDC hit
248  if (fabs(tdcTime) < timeWindowTDC) {
249  DTACHit *tacHitCombined = findMatch(tdcTime);
250  if (tacHitCombined == nullptr) {
251  // If there was no matching hit create a new one without FADC info
252  // push it into the list of hits for this factory.
253  tacHitCombined = new DTACHit();
254  tacHitCombined->setE(0.0);
255  tacHitCombined->setTimeFADC(
256  numeric_limits<double>::quiet_NaN());
257  tacHitCombined->setFADCPresent(false);
258  _data.push_back(tacHitCombined);
259  }
260 
261  tacHitCombined->setTDCPresent(true);
262  tacHitCombined->setT(tdcTime);
263 
264  // Check if the time-walk corrections are needed and adjust time time.
265  if (useTimeWalkCorrections && (tacHitCombined->getE() > 0.0)) {
266  // Correct for time walk
267  // The correction is the form t=t_tdc- C1 (A^C2 - A0^C2)
268  // double A = tacHitCombined->dE;
269  // double C1 = timewalk_parameters[id][1];
270  // double C2 = timewalk_parameters[id][2];
271  // double A0 = timewalk_parameters[id][3];
272  // T -= C1*(pow(A,C2) - pow(A0,C2));
273 
274  // Correct for timewalk using pulse peak instead of pulse integral
275  double A = tacHitCombined->getPulsePeak();
276 // double C0 = timeWalkParameters[0];
277  double C1 = timeWalkParameters[1];
278  double C2 = timeWalkParameters[2];
279  double A_THRESH = timeWalkParameters[3];
280  double A0 = timeWalkParameters[4];
281  // do correction
282  tdcTime -= C1
283  * (pow(A / A_THRESH, C2) - pow(A0 / A_THRESH, C2));
284  }
285  tacHitCombined->setT(tdcTime);
286 
287  tacHitCombined->AddAssociatedObject(tacTDCDigiHit);
288  }
289  }
290  return;
291 }
292 
293 //------------------
294 // Find a matching ADC hit and return pointer if found. It is assumed that no FADC hits have
295 // been inserted into the vector of the TAC hits, so looping through
296 // the vector will only take us through the ADC hits if anything at all.
297 // The ADC hit s are supposed to be taken care earlier but have NaN
298 // for TDC time information.
299 //------------------
301  DTACHit *bestMatchedHit = nullptr;
302 
303  // Loop over existing hits (from fADC) and look for a match
304  // in both the sector and the time.
305  for (auto& adcHit : _data) {
306  // only match to fADC hits, not bachelor TDC hits
307  if (!isfinite(adcHit->getTimeFADC()))
308  continue;
309 
310  double deltaT = fabs(adcHit->getTimeFADC() - tdcTime);
311  if (deltaT > timeDifferencInADCandTDC)
312  continue;
313 
314  if (bestMatchedHit != nullptr) {
315  if (deltaT < fabs(bestMatchedHit->getTimeFADC() - tdcTime))
316  bestMatchedHit = adcHit;
317  } else
318  bestMatchedHit = adcHit;
319  }
320  return bestMatchedHit;
321 }
322 
323 //------------------
324 // Reset_Data()
325 //------------------
327  // Clear _data vector
328  _data.clear();
329 }
330 
331 //------------------
332 // erun
333 //------------------
334 jerror_t DTACHit_factory::erun(void) {
335  return NOERROR;
336 }
337 
338 //------------------
339 // fini
340 //------------------
341 jerror_t DTACHit_factory::fini(void) {
342  return NOERROR;
343 }
344 
345 //std::string DTACHit_factory::SetTag(std::string tag) {
346 //
347 // char* newString = new char[tag.size()+1];
348 // std::copy( tag.begin(), tag.end(), newString );
349 // newString[tag.size()] = '\0';
350 // return string(tag_str);
351 //}
void setE(double e=0)
Definition: DTACHit.h:73
#define A0
Definition: src/md5.cpp:22
bool CheckFADC250_NoErrors(uint32_t QF) const
static set< int > announcedRuns
virtual jerror_t readCCDB(jana::JEventLoop *loop)
void setTimeTDC(double timeTdc=0)
Definition: DTACHit.h:121
double Convert_DigiTimeToNs_CAEN1290TDC(const JObject *locTDCDigiHit) const
static bool useTimeWalkCorrections
static bool checkErrorsOnFADC
double getPulsePeak() const
Definition: DTACHit.h:93
double getE() const
Definition: DTACHit.h:69
virtual void makeTDCHits(jana::JEventLoop *loop, uint64_t eventnumber)
void setTimeFADC(double timeFadc=0)
Definition: DTACHit.h:113
virtual jerror_t fini(void) override
Called after last event of last event source has been processed.
virtual void makeFADCHits(jana::JEventLoop *loop, uint64_t eventnumber)
virtual void Reset_Data(void)
virtual jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber) override
Called every event.
void setFADCPresent(bool fadcPresent=false)
Definition: DTACHit.h:129
void setPulsePeak(double pulsePeak=0)
Definition: DTACHit.h:97
virtual jerror_t erun(void) override
Called everytime run number changes, if brun has been called.
virtual DTACHit * findMatch(double tdcTime)
void setTDCPresent(bool tdcPresent=false)
Definition: DTACHit.h:137
static double timeWindowTDC
bool CheckFADC250_PedestalOK(uint32_t QF) const
virtual jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber) override
Called everytime a new run number is detected.
static double timeDifferencInADCandTDC
static TH1I * pedestal[nChan]
void setT(double t=0)
Definition: DTACHit.h:105
virtual jerror_t init(void) override
Called once at program start.
double getTimeFADC() const
Definition: DTACHit.h:109