Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HitRebuilderTAC.cc
Go to the documentation of this file.
1 /*
2  * HitRebuilderTAC.cc
3  *
4  * Created on: Jun 7, 2017
5  * Author: Hovanes Egiyan
6  */
7 
8 #include <cmath>
9 using namespace std;
10 
11 #include <TAC/HitRebuilderTAC.h>
12 
13 using namespace std;
14 
15 jerror_t HitRebuilderTAC::readCCDB( jana::JEventLoop* eventLoop) {
16  cout << "In HitRebuilderTAC::readCCDB() , reading calibration constants" << endl;
17  // load scale factors
18  map<string, double> scaleFactors;
19  if (eventLoop->GetCalib("/TAC/digi_scales", scaleFactors))
20  jout << "Error loading /TAC/digi_scales !" << endl;
21  // t_scale (TAC_ADC_SCALE)
22  if (scaleFactors.find("TAC_ADC_TSCALE") != scaleFactors.end())
23  timeScaleADC = adcTimeRescaleFactor * scaleFactors["TAC_ADC_TSCALE"];
24  else
25  jerr << "Unable to get TAC_ADC_TSCALE from /TAC/digi_scales !" << endl;
26 
27  // load base time offset
28  map<string, double> baseTimeOffsets;
29  // t_base (TAC_BASE_TIME_OFFSET)
30  if (eventLoop->GetCalib("/TAC/base_time_offset", baseTimeOffsets))
31  jout << "Error loading /TAC/base_time_offset !" << endl;
32  if (baseTimeOffsets.find("TAC_BASE_TIME_OFFSET") != baseTimeOffsets.end())
33  timeBaseADC = baseTimeOffsets["TAC_BASE_TIME_OFFSET"];
34  else
35  jerr
36  << "Unable to get TAC_BASE_TIME_OFFSET from /TAC/base_time_offset !"
37  << endl;
38 
39  // load constant tables
40  // adc_time_offsets (adc_timing_offsets)
41  if (eventLoop->GetCalib("/TAC/adc_timing_offsets", adcTimeOffset))
42  jout << "Error loading /TAC/adc_timing_offsets !" << endl;
43 
44  cout << "timeScaleADC is " << timeScaleADC << " , timeBaseADC is " << timeBaseADC <<
45  " , adcTimeOffset is " << adcTimeOffset << endl;
46 
47  return NOERROR;
48 }
49 
50 
52 //const vector<uint16_t>& DTACHit_Rebuild_factory::getSampleVector(
53 // const DTACHit* baseHit) {
54  set<const Df250WindowRawData*> rawDataSet;
55  set<const JObject*> alreadyChecked;
56  int maxDepth = 3;
57  baseHit->GetAssociatedAncestors(alreadyChecked, maxDepth, rawDataSet);
58 // cout << "Found " << associatedRawData.size() << " raw ancestors" << endl;
59 // for (auto& rawHit : associatedRawData) {
60 // cout << "There was a raw hit at " << hex << showbase << rawHit << dec
61 // << endl;
62 // }
63  // Make sure there is only one raw data object, otherwise return zero pointer
64  if (rawDataSet.size() != 1) {
65  stringstream errMsg;
66  errMsg << "DTACHit at " << showbase << baseHit
67  << " did not have Df250WindowRawData ancestors";
68  throw JException(errMsg.str());
69  }
70  return *rawDataSet.begin();
71 }
72 
73 double HitRebuilderTAC::getTimeFromRawData(const vector<uint16_t>& samples) {
74  if (samples.size() < 1)
75  return 0;
76  // Find the maximum by going through the raw data and comparing samples
77  pair<double, double> maxInfo(0, 0);
78  auto maxElementIter = std::max_element(samples.begin(), samples.end());
79  double peakLocation = std::distance(samples.begin(), maxElementIter);
80  // cout << "Maximum is at " << peakLocation << " bin" << endl;
81 
82  // CAlculated the weighted sum of the peak bin and the neighboring bins. This will
83  // be called the "hit time" and hopefully will consistent for the "fast pulses".
84  double weightedSum = 0;
85  double normalizationFactor = 0;
86 
87  auto minIter = maxElementIter - 1;
88  auto maxIter = maxElementIter + 3;
89  if (maxElementIter == samples.begin())
90  minIter = samples.begin();
91  if (maxElementIter >= (samples.end() - 2))
92  maxIter = samples.end();
93 
94  // cout << "minIter is at " << distance( samples.begin(), minIter )
95  // << " , maxIter is at " << distance( samples.begin(), maxIter )
96  // << endl;
97  for (auto& elemIter = minIter; elemIter != maxIter; elemIter++) {
98  double elemNumber = std::distance(samples.begin(), elemIter);
99  weightedSum += (elemNumber + 0.5) * (*elemIter);
100  normalizationFactor += (*elemIter);
101  }
102  // cout << "weighted sum is " << weightedSum << " norm factor is " << normalizationFactor << endl;
103  // If the normalization factor is non-zero, return the weighted average location
104  // of the tree bins around the highest bin.
105  if (normalizationFactor > 0)
106  peakLocation = weightedSum / normalizationFactor;
107 
108  // cout << "New peak location is " << peakLocation << endl;
109 
110  double pulseTime = timeScaleADC * peakLocation - adcTimeOffset
111  + timeBaseADC;
112 
113  // cout << "Rebuild pulse time will be " << pulseTime << endl;
114 
115  return pulseTime;
116 }
117 
119  vector<const DTACHit*>& baseHitVector) {
120  // Declare comparison functor for this TAC hits basedon pulse peak
121  static auto compareHits =
122  [](const DTACHit* lhs, const DTACHit* rhs ) ->
123  bool {return( (lhs!=nullptr)&(rhs!=nullptr ) ? fabs(lhs->getPulsePeak()) > fabs(rhs->getPulsePeak() ) : false );};
124 
125  // Sort the vector of the hits base on the peak height, with the heighest peak being in the first hit.
126  std::sort(baseHitVector.begin(), baseHitVector.end(), compareHits);
127 
128  rawDataPtrSet.clear();
129  vector<DTACHit*> newHitVector;
130  for (auto baseHit : baseHitVector) {
131  auto rawData = getRawData(baseHit);
132  // cout << "Raw data for " << baseHit << " is at " << rawData << endl;
133  // Check that no hit has already been made out of this raw data. This means that only
134  // one rebuild hit will be produced out of one window data. The fact that we sort the
135  // hits in the beginning should guarantee us that the useful hit is copied. This can be
136  // reprogrammed differently at a later time.
137  double newTime = baseHit->getT();
138  if (rawDataPtrSet.find(rawData) == rawDataPtrSet.end()) {
139  auto& rawDataVector = rawData->samples;
140  // auto rawDataVector = getSampleVector(baseHit);
141 
142  try {
143  newTime = getTimeFromRawData(rawDataVector);
144  } catch (const runtime_error& e) {
145  newTime = baseHit->getT();
146  newTime = -10;
147  }
148 // cout << "Received new time " << newTime << endl;
149 
150 // DTACHit *newHit = new DTACHit(*baseHit);
151  DTACHit* newHit = dynamic_cast<DTACHit*>(baseHit->Clone());
152 
153 // vector<pair<string,string>> outputStrings;
154 //
155 // baseHit->toStrings( outputStrings );
156 // cout << endl << "Old hit at " << baseHit << endl;
157 // for( auto& outputPair : outputStrings ) {
158 // cout << outputPair.first << " = " << outputPair.second << endl;
159 // }
160 // cout << endl << "New hit at " << newHit << endl;
161 // outputStrings.clear();
162 // newHit->toStrings( outputStrings );
163 // for( auto& outputPair : outputStrings ) {
164 // cout << outputPair.first << " = " << outputPair.second << endl;
165 // }
166 
167  newHit->setTimeFADC(newTime);
168  newHit->setT(newTime);
169 
170 // cout << "Times are : old = " << baseHit->getT() << " new = " << newHit->getT() << endl;
171 
172 // newHit->AddAssociatedObject(rawData);
173  newHit->AddAssociatedObject(baseHit);
174 
175  rawDataPtrSet.insert(rawData);
176 
177  newHitVector.push_back(newHit);
178  }
179  }
180  return newHitVector;
181 }
182 
jerror_t readCCDB(jana::JEventLoop *eventLoop)
virtual std::vector< DTACHit * > operator()(vector< const DTACHit * > &baseHitVector) override
double getPulsePeak() const
Definition: DTACHit.h:93
void setTimeFADC(double timeFadc=0)
Definition: DTACHit.h:113
TEllipse * e
virtual const Df250WindowRawData * getRawData(const DTACHit *baseHit) override
virtual double getTimeFromRawData(const std::vector< uint16_t > &samples) override
void setT(double t=0)
Definition: DTACHit.h:105