4 #include <JANA/JApplication.h>
17 if (enable_debug_output) {
18 bcal_points_tree =
new TTree(
"bcal_points_tree",
"");
19 bcal_points_tree->Branch(
"E",&E_tree,
"E/F");
20 bcal_points_tree->Branch(
"t_tdc",&t_tdc_tree,
"t_tdc/F");
21 bcal_points_tree->Branch(
"t_adc",&t_adc_tree,
"t_adc/F");
22 bcal_points_tree->Branch(
"t_tdc_corrected",&t_tdc_corrected_tree,
"t_tdc_corrected/F");
23 bcal_points_tree->Branch(
"t_adc_corrected",&t_adc_corrected_tree,
"t_adc_corrected/F");
24 bcal_points_tree->Branch(
"layer",&layer_tree,
"layer/I");
25 bcal_points_tree->Branch(
"end",&end_tree,
"end/O");
30 gPARMS->SetDefaultParameter(
"BCAL:USE_TDC", USE_TDC,
"Set to 1 to use TDC times");
33 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER;
34 static bool printMessage =
true;
35 pthread_mutex_lock(&print_mutex);
38 jout <<
"DBCALUnifiedHit_factory: Using TDC times when available." << endl;
41 jout <<
"DBCALUnifiedHit_factory: Using ADC times only." << endl;
45 pthread_mutex_unlock(&print_mutex);
53 vector<const DBCALGeometry *> BCALGeomVec;
54 eventLoop->Get(BCALGeomVec);
55 if(BCALGeomVec.size() == 0)
56 throw JException(
"Could not load DBCALGeometry object!");
57 dBCALGeom = BCALGeomVec[0];
60 JCalibration *jcalib = eventLoop->GetJCalibration();
62 vector<vector<float> > tdc_timewalk_table;
65 if(jcalib->Get(
"BCAL/timewalk_tdc_c4",tdc_timewalk_table)) {
66 jerr <<
"Error loading BCAL/timewalk_tdc !" << endl;
68 for (vector<vector<float> >::const_iterator iter = tdc_timewalk_table.begin();
69 iter != tdc_timewalk_table.end();
73 if (iter->size() != 9) {
74 cout <<
"DBCALUnifiedHit_factory: Wrong number of values in timewalk_tdc table (should be 9)" << endl;
78 int module = (int)((*iter)[0]+0.5);
79 int layer = (int)((*iter)[1]+0.5);
80 int sector = (int)((*iter)[2]+0.5);
81 int endi = (int)((*iter)[3]+0.5);
83 float c0 = (*iter)[4];
84 float c1 = (*iter)[5];
85 float c2 = (*iter)[6];
86 float c3 = (*iter)[7];
87 float thresh = (*iter)[8];
88 int cellId = dBCALGeom->cellId(module, layer, sector);
94 for (
int module=1; module<=dBCALGeom->GetBCAL_Nmodules(); module++) {
96 for (
int sector=1; sector<=4; sector++) {
98 int id = dBCALGeom->cellId(module,
layer, sector);
100 if (tdc_timewalk_map_c4.count(
readout_channel(
id,dBCALGeom->kUpstream)) != 1) {
101 cout <<
"DBCALUnifiedHit_factory: Channel missing in timewalk_tdc_table: "
102 << endl <<
" module " << module <<
" layer " <<
layer <<
" sector " << sector <<
" upstream" << endl;
105 if (tdc_timewalk_map_c4.count(
readout_channel(
id,dBCALGeom->kDownstream)) != 1) {
106 cout <<
"DBCALUnifiedHit_factory: Channel missing in timewalk_tdc_table: "
107 << endl <<
" module " << module <<
" layer " <<
layer <<
" sector " << sector <<
" downstream" << endl;
121 vector<const DBCALHit*> hits;
122 vector<const DBCALTDCHit*> tdc_hits;
125 if (hits.size() + tdc_hits.size() <= 0)
return NOERROR;
126 if (
VERBOSE>=3) jout << eventnumber <<
" " << hits.size() <<
" ADC hits, " << tdc_hits.size() <<
" TDC hits" << endl;
129 map<readout_channel, cellHits> cellHitMap;
130 for( vector<const DBCALHit*>::const_iterator hitPtr = hits.begin();
131 hitPtr != hits.end();
140 cellHitMap[
chan].hits.push_back(*hitPtr);
144 for(vector<const DBCALTDCHit*>::const_iterator hitPtr = tdc_hits.begin();
145 hitPtr != tdc_hits.end();
153 if( cellHitMap.find(chan) == cellHitMap.end() ){
157 cellHitMap[
chan].tdc_hits.push_back(*hitPtr);
161 for(map<readout_channel,cellHits>::const_iterator mapItr = cellHitMap.begin();
162 mapItr != cellHitMap.end();
167 int module = dBCALGeom->module(cellId);
168 int layer = dBCALGeom->layer(cellId);
169 int sector = dBCALGeom->sector(cellId);
171 const vector<const DBCALHit*> &hits = mapItr->second.hits;
172 const vector<const DBCALTDCHit*> &tdc_hits = mapItr->second.tdc_hits;
175 if (hits.size()==0) {
177 if (tdc_hits.size()!=0) {
178 static uint64_t Nwarnings = 0;
179 if(++Nwarnings <= 10) cout <<
"DBCALUnifiedHit_factory (event " << eventnumber <<
"): TDC hits without ADC hits" << endl;
180 if( Nwarnings == 10) cout <<
"DBCALUnifiedHit_factory: LAST WARNING (others will be suppressed)" <<endl;
187 unsigned int firstIndex = 0;
188 for(
unsigned int i=1; i<hits.size(); i++){
189 if (hits[i]->t < hits[firstIndex]->t) firstIndex = i;
190 if (
VERBOSE>=4 && hits.size()>=2) {
191 printf(
"DBCALUnifiedHit_factory event %5llu (%2i,%i,%i,%i) peak %4i, ADC %i %6.1f\n",
192 (
unsigned long long int)eventnumber, module, layer, sector, chan.
end, hits[i]->pulse_peak, i, hits[i]->t);
196 const DBCALHit* hit=hits[firstIndex];
197 float pulse_peak, E, t, t_ADC, t_TDC=0;
206 bool good_timewalk_params =
false;
207 if( tdc_timewalk_map_c4.find(chan) != tdc_timewalk_map_c4.end() ) {
210 tdc_coeff = tdc_timewalk_map_c4[
chan];
211 good_timewalk_params =
true;
214 bool haveTDChit =
false;
215 for (
unsigned int i=0; i<tdc_hits.size(); i++) {
218 float tdc_hit_t = tdc_hit->
t;
220 if(good_timewalk_params) {
222 float shifted_peak = pulse_peak+tdc_coeff.
c2;
223 if (shifted_peak>0) tdc_hit_t -= tdc_coeff.
c0 + tdc_coeff.
c1*pow(shifted_peak,tdc_coeff.
c3);
224 if (
VERBOSE>=4)
printf(
"tamewalk %f -> %f: (%f,%f,%f,%f)\n",tdc_hit->
t,tdc_hit_t,tdc_coeff.
c0,tdc_coeff.
c1,tdc_coeff.
c2,tdc_coeff.
c3);
227 float tdc_adc_diff = tdc_hit_t - t_ADC;
228 if (fabs(tdc_adc_diff) < fabs(t_diff)) {
234 printf(
"DBCALUnifiedHit_factory event %5llu (%2i,%i,%i,%i) peak %4.0f, ADC 0 %6.1f, TDC %i %6.1f diff %6.1f best %2i %6.1f\n",
235 (
unsigned long long int)eventnumber, module, layer, sector, chan.
end, pulse_peak, t_ADC, i, tdc_hit_t, tdc_adc_diff, goodTDCindex, t_diff);
238 if (
VERBOSE>=4 && !haveTDChit) {
241 printf(
"DBCALUnifiedHit_factory event %5llu (%2i,%i,%i,%i) peak %4.0f, ADC 0 %6.1f\n",
242 (
unsigned long long int)eventnumber, module, layer, sector, chan.
end, pulse_peak, hit->
t);
245 if (pulse_peak>tdc_coeff.
thresh && USE_TDC==1 && goodTDCindex>=0 && fabs(t_diff)<2 && good_timewalk_params) {
264 uhit->AddAssociatedObject(hit);
265 if (goodTDCindex>=0) uhit->AddAssociatedObject(tdc_hits[goodTDCindex]);
267 _data.push_back(uhit);
float t_TDC
Time of TDC hit that is closes t to the ADC time.
bool has_TDC_hit
Flag if the Unified Time is the TDC time.
static evioFileChannel * chan
float t_ADC
Time from fADC.
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
printf("string=%s", string)
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
float t
Unified time, obtained from ADC and/or TDC and used for further analysis.