File: | /home/sdobbs/work/clang/halld_recon/src/libraries/TRIGGER/DL1MCTrigger_factory.cc |
Warning: | line 617, column 2 Value stored to 'status' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
1 | // Test version V0 of the L1 trigger simulation, Aug 18, 2017 (asomov) |
2 | |
3 | #include <iostream> |
4 | #include <iomanip> |
5 | #include <cmath> |
6 | using namespace std; |
7 | |
8 | #include <JANA/JApplication.h> |
9 | #include <DAQ/DCODAROCInfo.h> |
10 | #include <DAQ/DL1Info.h> |
11 | #include <DANA/DStatusBits.h> |
12 | |
13 | #include <HDDM/DEventSourceHDDM.h> |
14 | |
15 | using namespace jana; |
16 | |
17 | #include "DL1MCTrigger_factory.h" |
18 | |
19 | #if HAVE_RCDB1 |
20 | #include "RCDB/Connection.h" |
21 | #include "RCDB/ConfigParser.h" |
22 | #endif |
23 | |
24 | static bool print_data_message = true; |
25 | |
26 | //------------------ |
27 | // init |
28 | //------------------ |
29 | jerror_t DL1MCTrigger_factory::init(void) |
30 | { |
31 | |
32 | debug = 0; |
33 | |
34 | if(debug){ |
35 | hfcal_gains = new TH1F("fcal_gains", "fcal_gains", 80, -1., 3.); |
36 | hfcal_gains2 = new TH2F("fcal_gains2","fcal_gains2", 71, -142., 142., 71, -142., 142.); |
37 | hfcal_ped = new TH1F("fcal_ped", "fcal_ped", 800, 0., 200.); |
38 | } |
39 | |
40 | BYPASS = 0; // default is to use trigger emulation |
41 | |
42 | // Default parameters for the main production trigger are taken from the |
43 | // spring run of 2017: 25 F + B > 45000 |
44 | |
45 | FCAL_ADC_PER_MEV = 3.73; |
46 | FCAL_CELL_THR = 65; |
47 | FCAL_EN_SC = 25; |
48 | FCAL_NSA = 10; |
49 | FCAL_NSB = 3; |
50 | FCAL_WINDOW = 10; |
51 | |
52 | BCAL_ADC_PER_MEV = 34.48276; // Not corrected energy |
53 | BCAL_CELL_THR = 20; |
54 | BCAL_EN_SC = 1; |
55 | BCAL_NSA = 19; |
56 | BCAL_NSB = 3; |
57 | BCAL_WINDOW = 20; |
58 | |
59 | FCAL_BCAL_EN = 45000; |
60 | |
61 | ST_ADC_PER_MEV = 1.; |
62 | ST_CELL_THR = 40; |
63 | ST_NSA = 10; |
64 | ST_NSB = 3; |
65 | ST_WINDOW = 10; |
66 | ST_NHIT = 1; |
67 | |
68 | BCAL_OFFSET = 2; |
69 | |
70 | SIMU_BASELINE = 1; |
71 | SIMU_GAIN = 1; |
72 | |
73 | |
74 | simu_baseline_fcal = 1; |
75 | simu_baseline_bcal = 1; |
76 | |
77 | simu_gain_fcal = 1; |
78 | simu_gain_bcal = 1; |
79 | |
80 | gPARMS->SetDefaultParameter("TRIG:BYPASS", BYPASS, |
81 | "Bypass trigger by hard coding physics bit"); |
82 | gPARMS->SetDefaultParameter("TRIG:FCAL_ADC_PER_MEV", FCAL_ADC_PER_MEV, |
83 | "FCAL energy calibration for the Trigger"); |
84 | gPARMS->SetDefaultParameter("TRIG:FCAL_CELL_THR", FCAL_CELL_THR, |
85 | "FCAL energy threshold per cell"); |
86 | gPARMS->SetDefaultParameter("TRIG:FCAL_EN_SC", FCAL_EN_SC, |
87 | "FCAL energy threshold"); |
88 | gPARMS->SetDefaultParameter("TRIG:FCAL_NSA", FCAL_NSA, |
89 | "FCAL NSA"); |
90 | gPARMS->SetDefaultParameter("TRIG:FCAL_NSB", FCAL_NSB, |
91 | "FCAL NSB"); |
92 | gPARMS->SetDefaultParameter("TRIG:FCAL_WINDOW", FCAL_WINDOW, |
93 | "FCAL GTP integration window"); |
94 | |
95 | gPARMS->SetDefaultParameter("TRIG:BCAL_ADC_PER_MEV", BCAL_ADC_PER_MEV, |
96 | "BCAL energy calibration for the Trigger"); |
97 | gPARMS->SetDefaultParameter("TRIG:BCAL_CELL_THR", BCAL_CELL_THR, |
98 | "BCAL energy threshold per cell"); |
99 | gPARMS->SetDefaultParameter("TRIG:BCAL_EN_SC", BCAL_EN_SC, |
100 | "BCAL energy threshold"); |
101 | gPARMS->SetDefaultParameter("TRIG:BCAL_NSA", BCAL_NSA, |
102 | "BCAL NSA"); |
103 | gPARMS->SetDefaultParameter("TRIG:BCAL_NSB", BCAL_NSB, |
104 | "BCAL NSB"); |
105 | gPARMS->SetDefaultParameter("TRIG:BCAL_WINDOW", BCAL_WINDOW, |
106 | "BCAL GTP integration window"); |
107 | |
108 | gPARMS->SetDefaultParameter("TRIG:ST_ADC_PER_MEV", ST_ADC_PER_MEV, |
109 | "ST energy calibration for the Trigger"); |
110 | gPARMS->SetDefaultParameter("TRIG:ST_CELL_THR", ST_CELL_THR, |
111 | "ST energy threshold per cell"); |
112 | gPARMS->SetDefaultParameter("TRIG:ST_NSA", ST_NSA, |
113 | "ST NSA"); |
114 | gPARMS->SetDefaultParameter("TRIG:ST_NSB", ST_NSB, |
115 | "ST NSB"); |
116 | gPARMS->SetDefaultParameter("TRIG:ST_WINDOW", ST_WINDOW, |
117 | "ST window for merging hits (GTP)"); |
118 | gPARMS->SetDefaultParameter("TRIG:ST_NHIT", ST_NHIT, |
119 | "Number of hits in ST"); |
120 | |
121 | gPARMS->SetDefaultParameter("TRIG:FCAL_BCAL_EN", FCAL_BCAL_EN, |
122 | "Energy threshold for the FCAL & BCAL trigger"); |
123 | |
124 | gPARMS->SetDefaultParameter("TRIG:BCAL_OFFSET", BCAL_OFFSET, |
125 | "Timing offset between BCAL and FCAL energies at GTP (sampels)"); |
126 | |
127 | |
128 | // Allows to switch off gain and baseline fluctuations |
129 | gPARMS->SetDefaultParameter("TRIG:SIMU_BASELINE", SIMU_BASELINE, |
130 | "Enable simulation of pedestal variations"); |
131 | |
132 | gPARMS->SetDefaultParameter("TRIG:SIMU_GAIN", SIMU_GAIN, |
133 | "Enable simulation of gain variations"); |
134 | |
135 | |
136 | BCAL_ADC_PER_MEV_CORRECT = 22.7273; |
137 | |
138 | pedestal_sigma = 1.2; |
139 | |
140 | time_shift = 100; |
141 | |
142 | time_min = 0; |
143 | time_max = (sample - 1)*max_adc_bins; |
144 | |
145 | vector< vector<double > > fcal_gains_temp(DFCALGeometry::kBlocksTall, |
146 | vector<double>(DFCALGeometry::kBlocksWide)); |
147 | vector< vector<double > > fcal_pedestals_temp(DFCALGeometry::kBlocksTall, |
148 | vector<double>(DFCALGeometry::kBlocksWide)); |
149 | |
150 | fcal_gains = fcal_gains_temp; |
151 | fcal_pedestals = fcal_pedestals_temp; |
152 | |
153 | return NOERROR; |
154 | } |
155 | |
156 | |
157 | //------------------ |
158 | // brun |
159 | //------------------ |
160 | jerror_t DL1MCTrigger_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) |
161 | { |
162 | |
163 | int use_rcdb = 1; |
164 | |
165 | int status = 0; |
166 | |
167 | fcal_trig_mask.clear(); |
168 | bcal_trig_mask.clear(); |
169 | |
170 | triggers_enabled.clear(); |
171 | |
172 | // Only print messages for one thread whenever run number change |
173 | |
174 | static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
175 | static set<int> runs_announced; |
176 | pthread_mutex_lock(&print_mutex); |
177 | bool print_messages = false; |
178 | |
179 | if(runs_announced.find(runnumber) == runs_announced.end()){ |
180 | print_messages = true; |
181 | runs_announced.insert(runnumber); |
182 | } |
183 | |
184 | pthread_mutex_unlock(&print_mutex); |
185 | |
186 | // Don't use RCDB for mc_generic: ideal MC simulation |
187 | |
188 | string JANA_CALIB_CONTEXT = ""; |
189 | |
190 | if(getenv("JANA_CALIB_CONTEXT") != NULL__null ){ |
191 | JANA_CALIB_CONTEXT = getenv("JANA_CALIB_CONTEXT"); |
192 | if(JANA_CALIB_CONTEXT.find("mc_generic") != string::npos){ |
193 | use_rcdb = 0; |
194 | // Don't simulate baseline fluctuations for mc_generic |
195 | simu_baseline_fcal = 0; |
196 | simu_baseline_bcal = 0; |
197 | // Don't simulate gain fluctuations for mc_generic |
198 | simu_gain_fcal = 0; |
199 | simu_gain_bcal = 0; |
200 | } |
201 | } |
202 | |
203 | // runnumber = 30942; |
204 | |
205 | if(use_rcdb == 1){ |
206 | status = Read_RCDB(runnumber, print_messages); |
207 | if(print_messages) PrintTriggers(); |
208 | } |
209 | |
210 | |
211 | if( (use_rcdb == 0) || (status > 0) || (triggers_enabled.size() == 0)){ |
212 | |
213 | // Simulate FCAL & BCAL main production trigger only |
214 | |
215 | trigger_conf trig_tmp; |
216 | trig_tmp.bit = 0; |
217 | trig_tmp.gtp.fcal = FCAL_EN_SC; |
218 | trig_tmp.gtp.bcal = BCAL_EN_SC; |
219 | trig_tmp.gtp.en_thr = FCAL_BCAL_EN; |
220 | trig_tmp.gtp.fcal_min = 200; |
221 | triggers_enabled.push_back(trig_tmp); |
222 | |
223 | if(print_messages) |
224 | cout << " Do not use RCDB for the trigger simulation. Default (spring 2017) trigger settings are used " << endl; |
225 | } |
226 | |
227 | |
228 | // extract the FCAL Geometry |
229 | vector<const DFCALGeometry*> fcalGeomVect; |
230 | eventLoop->Get( fcalGeomVect ); |
231 | if (fcalGeomVect.size() < 1) |
232 | return OBJECT_NOT_AVAILABLE; |
233 | const DFCALGeometry& fcalGeom = *(fcalGeomVect[0]); |
234 | |
235 | if(print_messages) jout << "In DL1MCTrigger_factory, loading constants..." << endl; |
236 | |
237 | vector< double > fcal_gains_ch; |
238 | vector< double > fcal_pedestals_ch; |
239 | |
240 | if (eventLoop->GetCalib("/FCAL/gains", fcal_gains_ch)){ |
241 | jout << "DL1MCTrigger_factory: Error loading /FCAL/gains !" << endl; |
242 | // Load default values of gains if CCDB table is not found |
243 | for(int ii = 0; ii < DFCALGeometry::kBlocksTall; ii++){ |
244 | for(int jj = 0; jj < DFCALGeometry::kBlocksWide; jj++){ |
245 | fcal_gains[ii][jj] = 1.; |
246 | } |
247 | } |
248 | } else { |
249 | LoadFCALConst(fcal_gains, fcal_gains_ch, fcalGeom); |
250 | |
251 | if(debug){ |
252 | for(int ch = 0; ch < (int)fcal_gains_ch.size(); ch++){ |
253 | int row = fcalGeom.row(ch); |
254 | int col = fcalGeom.column(ch); |
255 | // Sanity check for regular FCAL (row,col) ranges (anticipating |
256 | // future upgrade to FCAL to include insert) |
257 | if(fcalGeom.isBlockActive(row,col)&&row<DFCALGeometry::kBlocksTall |
258 | && col<DFCALGeometry::kBlocksWide){ |
259 | hfcal_gains->Fill(fcal_gains[row][col]); |
260 | DVector2 aaa = fcalGeom.positionOnFace(row,col); |
261 | hfcal_gains2->Fill(float(aaa.X()), float(aaa.Y()), fcal_gains[row][col]); |
262 | // cout << aaa.X() << " " << aaa.Y() << endl; |
263 | |
264 | } |
265 | } |
266 | } |
267 | |
268 | } |
269 | |
270 | if (eventLoop->GetCalib("/FCAL/pedestals", fcal_pedestals_ch)){ |
271 | jout << "DL1MCTrigger_factory: Error loading /FCAL/pedestals !" << endl; |
272 | // Load default values of pedestals if CCDB table is not found |
273 | for(int ii = 0; ii < DFCALGeometry::kBlocksTall; ii++){ |
274 | for(int jj = 0; jj < DFCALGeometry::kBlocksWide; jj++){ |
275 | fcal_pedestals[ii][jj] = 100.; |
276 | } |
277 | } |
278 | } else { |
279 | LoadFCALConst(fcal_pedestals, fcal_pedestals_ch, fcalGeom); |
280 | |
281 | if(debug){ |
282 | for(int ch = 0; ch < (int)fcal_gains_ch.size(); ch++){ |
283 | int row = fcalGeom.row(ch); |
284 | int col = fcalGeom.column(ch); |
285 | // Sanity check for regular FCAL (row,col) ranges (anticipating |
286 | // future upgrade to FCAL to include insert) |
287 | if(fcalGeom.isBlockActive(row,col)&&row<DFCALGeometry::kBlocksTall |
288 | && col<DFCALGeometry::kBlocksWide){ |
289 | hfcal_ped->Fill(fcal_pedestals[row][col]); |
290 | } |
291 | } |
292 | } |
293 | |
294 | } |
295 | |
296 | if(!SIMU_BASELINE){ |
297 | simu_baseline_fcal = 0; |
298 | simu_baseline_bcal = 0; |
299 | } |
300 | |
301 | if(!SIMU_GAIN){ |
302 | simu_gain_fcal = 0; |
303 | simu_gain_bcal = 0; |
304 | } |
305 | |
306 | if(debug){ |
307 | for(int ii = 0; ii < 100; ii++){ |
308 | cout << " Channel = " << ii << " Value = " << |
309 | fcal_gains_ch[ii] << endl; |
310 | } |
311 | } |
312 | |
313 | |
314 | return NOERROR; |
315 | } |
316 | |
317 | //------------------ |
318 | // evnt |
319 | //------------------ |
320 | jerror_t DL1MCTrigger_factory::evnt(JEventLoop *loop, uint64_t eventnumber){ |
321 | |
322 | if(BYPASS) { |
323 | DL1MCTrigger *trigger = new DL1MCTrigger; |
324 | trigger->trig_mask = 1; |
325 | _data.push_back(trigger); |
326 | return NOERROR; |
327 | } |
328 | |
329 | int l1_found = 1; |
330 | |
331 | int status = 0; |
332 | |
333 | fcal_signal_hits.clear(); |
334 | bcal_signal_hits.clear(); |
335 | |
336 | fcal_merged_hits.clear(); |
337 | bcal_merged_hits.clear(); |
338 | |
339 | memset(fcal_ssp,0,sizeof(fcal_ssp)); |
340 | memset(fcal_gtp,0,sizeof(fcal_gtp)); |
341 | |
342 | memset(bcal_ssp,0,sizeof(bcal_ssp)); |
343 | memset(bcal_gtp,0,sizeof(bcal_gtp)); |
344 | |
345 | |
346 | vector<const DFCALHit*> fcal_hits; |
347 | vector<const DBCALHit*> bcal_hits; |
348 | |
349 | loop->Get(fcal_hits); |
350 | loop->Get(bcal_hits); |
351 | |
352 | DRandom2 gDRandom(0); // declared extern in DRandom2.h |
353 | |
354 | // This is temporary, to allow this simulation to be run on data |
355 | // to help out with trigger efficiency studies - sdobbs (Aug. 26, 2020) |
356 | if( loop->GetJEvent().GetStatusBit(kSTATUS_EVIO) ){ |
357 | if(print_data_message) { |
358 | jout << "WARNING: Running L1 trigger simulation on EVIO data" << endl; |
359 | print_data_message = false; |
360 | } |
361 | |
362 | // for data, don't add in baseline shifts, since they already exist |
363 | simu_baseline_fcal = 0; |
364 | simu_baseline_bcal = 0; |
365 | } else { |
366 | // Initialize random number generator |
367 | // Read seeds from hddm file |
368 | // Generate seeds according to the event number if they are not stored in hddm |
369 | // The proceedure is consistent with the mcsmear |
370 | |
371 | UInt_t seed1 = 0; |
372 | UInt_t seed2 = 0; |
373 | UInt_t seed3 = 0; |
374 | |
375 | GetSeeds(loop, eventnumber, seed1, seed2, seed3); |
376 | |
377 | gDRandom.SetSeeds(seed1, seed2, seed3); |
378 | } |
379 | |
380 | // cout << endl; |
381 | // cout << " Event = " << eventnumber << endl; |
382 | // cout << " Seed 1: " << seed1 << endl; |
383 | // cout << " Seed 2: " << seed2 << endl; |
384 | // cout << " Seed 3: " << seed3 << endl; |
385 | // cout << endl; |
386 | |
387 | |
388 | DL1MCTrigger *trigger = new DL1MCTrigger; |
389 | |
390 | // FCAL energy sum |
391 | double fcal_hit_en = 0; |
392 | |
393 | for (unsigned int ii = 0; ii < fcal_hits.size(); ii++){ |
394 | int row = fcal_hits[ii]->row; |
395 | int col = fcal_hits[ii]->column; |
396 | |
397 | // Shift time to simulate pile up hits |
398 | double time = fcal_hits[ii]->t + time_shift; |
399 | if((time < time_min) || (time > time_max)){ |
400 | continue; |
401 | } |
402 | |
403 | // Check channels masked for trigger |
404 | int ch_masked = 0; |
405 | for(unsigned int jj = 0; jj < fcal_trig_mask.size(); jj++){ |
406 | if( (row == fcal_trig_mask[jj].row) && (col == fcal_trig_mask[jj].col)){ |
407 | ch_masked = 1; |
408 | break; |
409 | } |
410 | } |
411 | |
412 | if(ch_masked == 0){ |
413 | fcal_hit_en += fcal_hits[ii]->E; |
414 | |
415 | fcal_signal fcal_tmp; |
416 | |
417 | fcal_tmp.row = row; |
418 | fcal_tmp.column = col; |
419 | |
420 | fcal_tmp.energy = fcal_hits[ii]->E; |
421 | fcal_tmp.time = time; |
422 | memset(fcal_tmp.adc_amp,0,sizeof(fcal_tmp.adc_amp)); |
423 | memset(fcal_tmp.adc_en, 0,sizeof(fcal_tmp.adc_en)); |
424 | |
425 | double fcal_adc_en = fcal_tmp.energy*FCAL_ADC_PER_MEV*1000; |
426 | |
427 | // Account for gain fluctuations |
428 | if(simu_gain_fcal && row<DFCALGeometry::kBlocksTall |
429 | && col<DFCALGeometry::kBlocksWide){ |
430 | |
431 | double gain = fcal_gains[row][col]; |
432 | |
433 | fcal_adc_en *= gain; |
434 | } |
435 | |
436 | status = SignalPulse(fcal_adc_en, fcal_tmp.time, fcal_tmp.adc_en, 1); |
437 | status = 0; |
438 | |
439 | fcal_signal_hits.push_back(fcal_tmp); |
440 | } |
441 | |
442 | } |
443 | |
444 | |
445 | // Merge FCAL hits |
446 | for(unsigned int ii = 0; ii < fcal_signal_hits.size(); ii++){ |
447 | |
448 | if(fcal_signal_hits[ii].merged == 1) continue; |
449 | |
450 | fcal_signal fcal_tmp; |
451 | fcal_tmp.row = fcal_signal_hits[ii].row; |
452 | fcal_tmp.column = fcal_signal_hits[ii].column; |
453 | fcal_tmp.energy = 0.; |
454 | fcal_tmp.time = 0.; |
455 | for(int kk = 0; kk < sample; kk++) |
456 | fcal_tmp.adc_en[kk] = fcal_signal_hits[ii].adc_en[kk]; |
457 | |
458 | for(unsigned int jj = ii + 1; jj < fcal_signal_hits.size(); jj++){ |
459 | if((fcal_signal_hits[ii].row == fcal_signal_hits[jj].row) && |
460 | (fcal_signal_hits[ii].column == fcal_signal_hits[jj].column)){ |
461 | |
462 | fcal_signal_hits[jj].merged = 1; |
463 | |
464 | for(int kk = 0; kk < sample; kk++) |
465 | fcal_tmp.adc_en[kk] += fcal_signal_hits[jj].adc_en[kk]; |
466 | } |
467 | } |
468 | |
469 | fcal_merged_hits.push_back(fcal_tmp); |
470 | } |
471 | |
472 | // Add baseline fluctuations for channels with hits |
473 | if(simu_baseline_fcal){ |
474 | for(unsigned int ii = 0; ii < fcal_merged_hits.size(); ii++){ |
475 | int row = fcal_merged_hits[ii].row; |
476 | int column = fcal_merged_hits[ii].column; |
477 | double pedestal = 100.0; |
478 | if (row<DFCALGeometry::kBlocksTall |
479 | && column<DFCALGeometry::kBlocksWide){ |
480 | pedestal=fcal_pedestals[row][column]; |
481 | } |
482 | AddBaseline(fcal_merged_hits[ii].adc_en, pedestal, gDRandom); |
483 | } |
484 | } |
485 | |
486 | |
487 | // Digitize |
488 | for(unsigned int ii = 0; ii < fcal_merged_hits.size(); ii++){ |
489 | Digitize(fcal_merged_hits[ii].adc_en,fcal_merged_hits[ii].adc_amp); |
490 | // cout << " Digitize " << fcal_merged_hits[ii].adc_en[sample - 3] |
491 | // << " " << fcal_merged_hits[ii].adc_amp[sample - 3] << endl; |
492 | } |
493 | |
494 | |
495 | int fcal_hit_adc_en = 0; |
496 | |
497 | for(unsigned int ii = 0; ii < fcal_merged_hits.size(); ii++) |
498 | for(int jj = 0; jj < sample; jj++) |
499 | fcal_hit_adc_en += fcal_merged_hits[ii].adc_amp[jj]; |
500 | |
501 | |
502 | status += FADC_SSP(fcal_merged_hits, 1); |
503 | |
504 | status += GTP(1); |
505 | |
506 | |
507 | |
508 | // BCAL energy sum |
509 | double bcal_hit_en = 0; |
510 | |
511 | for (unsigned int ii = 0; ii < bcal_hits.size(); ii++){ |
512 | |
513 | // Shift time to simulate pile up hits |
514 | double time = bcal_hits[ii]->t + time_shift; |
515 | if((time < time_min) || (time > time_max)){ |
516 | continue; |
517 | } |
518 | |
519 | int module = bcal_hits[ii]->module; |
520 | int layer = bcal_hits[ii]->layer; |
521 | int sector = bcal_hits[ii]->sector; |
522 | int end = bcal_hits[ii]->end; |
523 | |
524 | // Check trigger masks |
525 | int ch_masked = 0; |
526 | |
527 | for(unsigned int jj = 0; jj < bcal_trig_mask.size(); jj++){ |
528 | if( (module == bcal_trig_mask[jj].module) && (layer == bcal_trig_mask[jj].layer) && |
529 | (sector == bcal_trig_mask[jj].sector) && (end == bcal_trig_mask[jj].end) ){ |
530 | ch_masked = 1; |
531 | break; |
532 | } |
533 | } |
534 | |
535 | |
536 | if(ch_masked == 0){ |
537 | bcal_hit_en += bcal_hits[ii]->E; |
538 | |
539 | bcal_signal bcal_tmp; |
540 | bcal_tmp.module = module; |
541 | bcal_tmp.layer = layer; |
542 | bcal_tmp.sector = sector; |
543 | bcal_tmp.end = end; |
544 | |
545 | bcal_tmp.energy = bcal_hits[ii]->E; |
546 | bcal_tmp.time = time; |
547 | memset(bcal_tmp.adc_amp,0,sizeof(bcal_tmp.adc_amp)); |
548 | memset(bcal_tmp.adc_en, 0,sizeof(bcal_tmp.adc_en)); |
549 | |
550 | double bcal_adc_en = bcal_tmp.energy*BCAL_ADC_PER_MEV*1000; |
551 | |
552 | status = SignalPulse(bcal_adc_en, bcal_tmp.time, bcal_tmp.adc_en, 2); |
553 | status = 0; |
554 | |
555 | bcal_signal_hits.push_back(bcal_tmp); |
556 | } |
557 | } |
558 | |
559 | |
560 | // Merge BCAL hits |
561 | for(unsigned int ii = 0; ii < bcal_signal_hits.size(); ii++){ |
562 | |
563 | if(bcal_signal_hits[ii].merged == 1) continue; |
564 | |
565 | bcal_signal bcal_tmp; |
566 | bcal_tmp.module = bcal_signal_hits[ii].module; |
567 | bcal_tmp.layer = bcal_signal_hits[ii].layer; |
568 | bcal_tmp.sector = bcal_signal_hits[ii].sector; |
569 | bcal_tmp.end = bcal_signal_hits[ii].end; |
570 | |
571 | bcal_tmp.energy = 0.; |
572 | bcal_tmp.time = 0.; |
573 | |
574 | for(int kk = 0; kk < sample; kk++) |
575 | bcal_tmp.adc_en[kk] = bcal_signal_hits[ii].adc_en[kk]; |
576 | |
577 | for(unsigned int jj = ii + 1; jj < bcal_signal_hits.size(); jj++){ |
578 | if((bcal_signal_hits[ii].module == bcal_signal_hits[jj].module) && |
579 | (bcal_signal_hits[ii].layer == bcal_signal_hits[jj].layer) && |
580 | (bcal_signal_hits[ii].sector == bcal_signal_hits[jj].sector) && |
581 | (bcal_signal_hits[ii].end == bcal_signal_hits[jj].end)){ |
582 | |
583 | bcal_signal_hits[jj].merged = 1; |
584 | |
585 | for(int kk = 0; kk < sample; kk++) |
586 | bcal_tmp.adc_en[kk] += bcal_signal_hits[jj].adc_en[kk]; |
587 | } |
588 | } |
589 | |
590 | bcal_merged_hits.push_back(bcal_tmp); |
591 | } |
592 | |
593 | |
594 | // Add baseline fluctuations for channels with hits |
595 | if(simu_baseline_bcal){ |
596 | for(unsigned int ii = 0; ii < bcal_merged_hits.size(); ii++){ |
597 | // Assume that all BCAL pedestals are 100 |
598 | double pedestal = TRIG_BASELINE; |
599 | AddBaseline(bcal_merged_hits[ii].adc_en, pedestal, gDRandom); |
600 | } |
601 | } |
602 | |
603 | |
604 | // Digitize |
605 | for(unsigned int ii = 0; ii < bcal_merged_hits.size(); ii++) |
606 | Digitize(bcal_merged_hits[ii].adc_en,bcal_merged_hits[ii].adc_amp); |
607 | |
608 | |
609 | int bcal_hit_adc_en = 0; |
610 | for(unsigned int ii = 0; ii < bcal_merged_hits.size(); ii++) |
611 | for(int jj = 0; jj < sample; jj++) |
612 | bcal_hit_adc_en += bcal_merged_hits[ii].adc_amp[jj]; |
613 | |
614 | |
615 | status = FADC_SSP(bcal_merged_hits, 2); |
616 | |
617 | status = GTP(2); |
Value stored to 'status' is never read | |
618 | |
619 | // Search for triggers |
620 | |
621 | l1_found = FindTriggers(trigger); |
622 | |
623 | if(l1_found){ |
624 | |
625 | int fcal_gtp_max = 0; |
626 | int bcal_gtp_max = 0; |
627 | |
628 | for(unsigned int ii = 0; ii < sample; ii++){ |
629 | if(fcal_gtp[ii] > fcal_gtp_max) fcal_gtp_max = fcal_gtp[ii]; |
630 | if(bcal_gtp[ii] > bcal_gtp_max) bcal_gtp_max = bcal_gtp[ii]; |
631 | } |
632 | |
633 | trigger->fcal_en = fcal_hit_en; |
634 | trigger->fcal_adc = fcal_hit_adc_en; |
635 | trigger->fcal_adc_en = fcal_hit_adc_en/FCAL_ADC_PER_MEV/1000.; |
636 | trigger->fcal_gtp = fcal_gtp_max; |
637 | trigger->fcal_gtp_en = fcal_gtp_max/FCAL_ADC_PER_MEV/1000.; |
638 | |
639 | trigger->bcal_en = bcal_hit_en; |
640 | trigger->bcal_adc = bcal_hit_adc_en; |
641 | trigger->bcal_adc_en = bcal_hit_adc_en/BCAL_ADC_PER_MEV_CORRECT/2./1000.; |
642 | trigger->bcal_gtp = bcal_gtp_max; |
643 | trigger->bcal_gtp_en = bcal_gtp_max/BCAL_ADC_PER_MEV_CORRECT/2./1000.; |
644 | |
645 | _data.push_back(trigger); |
646 | |
647 | } else{ |
648 | delete trigger; |
649 | } |
650 | |
651 | return NOERROR; |
652 | } |
653 | |
654 | //------------------ |
655 | // erun |
656 | //------------------ |
657 | jerror_t DL1MCTrigger_factory::erun(void) |
658 | { |
659 | return NOERROR; |
660 | } |
661 | |
662 | //------------------ |
663 | // fini |
664 | //------------------ |
665 | jerror_t DL1MCTrigger_factory::fini(void) |
666 | { |
667 | return NOERROR; |
668 | } |
669 | |
670 | //********************* |
671 | // Read RCDB |
672 | //********************* |
673 | |
674 | int DL1MCTrigger_factory::Read_RCDB(int32_t runnumber, bool print_messages) |
675 | { |
676 | |
677 | #if HAVE_RCDB1 |
678 | |
679 | vector<const DTranslationTable*> ttab; |
680 | eventLoop->Get(ttab); |
681 | |
682 | vector<string> SectionNames = {"TRIGGER", "GLOBAL", "FCAL", "BCAL", "TOF", "ST", "TAGH", |
683 | "TAGM", "PS", "PSC", "TPOL", "CDC", "FDC"}; |
684 | |
685 | string RCDB_CONNECTION; |
686 | |
687 | if( getenv("RCDB_CONNECTION")!= NULL__null ) |
688 | RCDB_CONNECTION = getenv("RCDB_CONNECTION"); |
689 | else |
690 | RCDB_CONNECTION = "mysql://rcdb@hallddb.jlab.org/rcdb"; // default to outward-facing MySQL DB |
691 | |
692 | |
693 | rcdb::Connection connection(RCDB_CONNECTION); |
694 | |
695 | auto rtvsCnd = connection.GetCondition(runnumber, "rtvs"); |
696 | |
697 | if( !rtvsCnd ) { |
698 | cout<<"'rtvs' condition is not set for run " << runnumber << endl; |
699 | return 2; |
700 | } |
701 | |
702 | |
703 | auto json = rtvsCnd->ToJsonDocument(); // The CODA rtvs is serialized as JSon dictionary. |
704 | string fileName(json["%(config)"].GetString()); // We need item with name '%(config)' |
705 | |
706 | auto file = connection.GetFile(runnumber, fileName); |
707 | |
708 | if(!file) { // If there is no such file, null is returned |
709 | cout<<"File with name: "<< fileName |
710 | <<" doesn't exist (not associated) with run: "<< runnumber << endl; |
711 | return 3; |
712 | } |
713 | |
714 | string fileContent = file->GetContent(); // Get file content |
715 | auto result = rcdb::ConfigParser::Parse(fileContent, SectionNames); // Parse it! |
716 | |
717 | |
718 | // Pulse parameters for trigger |
719 | |
720 | auto trig_thr = result.Sections["FCAL"].NameValues["FADC250_TRIG_THR"]; |
721 | auto trig_nsb = result.Sections["FCAL"].NameValues["FADC250_TRIG_NSB"]; |
722 | auto trig_nsa = result.Sections["FCAL"].NameValues["FADC250_TRIG_NSA"]; |
723 | |
724 | if(trig_thr.size() > 0){ |
725 | FCAL_CELL_THR = stoi(trig_thr); |
726 | if(FCAL_CELL_THR < 0) FCAL_CELL_THR = 0; |
727 | } |
728 | |
729 | if(trig_nsb.size() > 0) |
730 | FCAL_NSB = stoi(trig_nsb); |
731 | |
732 | if(trig_thr.size() > 0) |
733 | FCAL_NSA = stoi(trig_nsa); |
734 | |
735 | trig_thr = result.Sections["BCAL"].NameValues["FADC250_TRIG_THR"]; |
736 | trig_nsb = result.Sections["BCAL"].NameValues["FADC250_TRIG_NSB"]; |
737 | trig_nsa = result.Sections["BCAL"].NameValues["FADC250_TRIG_NSA"]; |
738 | |
739 | if(trig_thr.size() > 0){ |
740 | BCAL_CELL_THR = stoi(trig_thr); |
741 | if(BCAL_CELL_THR < 0) BCAL_CELL_THR = 0; |
742 | } |
743 | |
744 | if(trig_nsb.size() > 0) |
745 | BCAL_NSB = stoi(trig_nsb); |
746 | |
747 | if(trig_thr.size() > 0) |
748 | BCAL_NSA = stoi(trig_nsa); |
749 | |
750 | |
751 | // List of enabled GTP equations and triggers |
752 | vector<vector<string>> triggerTypes; |
753 | |
754 | for(auto row : result.Sections["TRIGGER"].Rows) { |
755 | |
756 | if(row[0] == "TRIG_TYPE") { |
757 | |
758 | if(row.size() >= 9){ |
759 | triggerTypes.push_back(row); // The line in config file starts with TRIG_TYPE |
760 | } else { |
761 | cout << " Cannot parse TRIG_TYPE. Insufficient number of parameters " << row.size() << endl; |
762 | } |
763 | } |
764 | |
765 | |
766 | // Integration windows for BCAL and FCAL |
767 | if(row[0] == "TRIG_EQ") { |
768 | if(row.size() >= 5){ |
769 | if(stoi(row[4]) == 1){ // Trigger enabled |
770 | if(row[1] == "FCAL") |
771 | FCAL_WINDOW = stoi(row[3]); |
772 | if(row[1] == "BCAL_E") |
773 | BCAL_WINDOW = stoi(row[3]); |
774 | } |
775 | } |
776 | } |
777 | } |
778 | |
779 | |
780 | |
781 | |
782 | for(int ii = 0; ii < 32; ii++){ |
783 | |
784 | int trig_found = 0; |
785 | |
786 | trigger_conf trigger_tmp; |
787 | |
788 | memset(&trigger_tmp,0,sizeof(trigger_tmp)); |
789 | |
790 | trigger_tmp.bit = ii; |
791 | |
792 | for(unsigned int jj = 0; jj < triggerTypes.size(); jj++){ |
793 | |
794 | if(triggerTypes[jj].size() < 9) continue; // Minimum 9 parameters are required in the config file |
795 | |
796 | int bit = stoi(triggerTypes[jj][8]); // Trigger bit |
797 | |
798 | if( bit == ii){ |
799 | |
800 | // FCAL, BCAL triggers |
801 | if(triggerTypes[jj][1] == "BFCAL"){ |
802 | |
803 | int fcal = 0; |
804 | int bcal = 0; |
805 | |
806 | if(triggerTypes[jj][4].size() > 0) fcal = stoi(triggerTypes[jj][4]); |
807 | if(triggerTypes[jj][5].size() > 0) bcal = stoi(triggerTypes[jj][5]); |
808 | |
809 | |
810 | if( (fcal > 0) && (bcal > 0)){ // FCAL & BCAL |
811 | trigger_tmp.type = 0x3; |
812 | trigger_tmp.gtp.fcal = fcal; |
813 | trigger_tmp.gtp.bcal = bcal; |
814 | if(triggerTypes[jj][6].size() > 0) trigger_tmp.gtp.en_thr = stoi(triggerTypes[jj][6]); |
815 | |
816 | if(triggerTypes[jj].size() > 9){ |
817 | if((triggerTypes[jj].size() >= 10) && (triggerTypes[jj][9].size() > 0)) trigger_tmp.gtp.fcal_min = stoi(triggerTypes[jj][9]); |
818 | if((triggerTypes[jj].size() >= 11) && (triggerTypes[jj][10].size() > 0)) trigger_tmp.gtp.fcal_max = stoi(triggerTypes[jj][10]); |
819 | if((triggerTypes[jj].size() >= 12) && (triggerTypes[jj][11].size() > 0)) trigger_tmp.gtp.bcal_min = stoi(triggerTypes[jj][11]); |
820 | if((triggerTypes[jj].size() >= 13) && (triggerTypes[jj][12].size() > 0)) trigger_tmp.gtp.bcal_max = stoi(triggerTypes[jj][12]); |
821 | } |
822 | |
823 | trig_found++; |
824 | } else if ((fcal > 0) && (bcal == 0)){ // FCAL only |
825 | trigger_tmp.type = 0x1; |
826 | trigger_tmp.gtp.fcal = fcal; |
827 | if(triggerTypes[jj][6].size() > 0) trigger_tmp.gtp.en_thr = stoi(triggerTypes[jj][6]); |
828 | if(triggerTypes[jj].size() > 9){ |
829 | if(triggerTypes[jj][9].size() > 0) trigger_tmp.gtp.fcal_min = stoi(triggerTypes[jj][9]); |
830 | if(triggerTypes[jj][10].size() > 0) trigger_tmp.gtp.fcal_max = stoi(triggerTypes[jj][10]); |
831 | } |
832 | trig_found++; |
833 | } else if ((bcal > 0) && (fcal == 0)){ // BCAL only |
834 | trigger_tmp.type = 0x2; |
835 | trigger_tmp.gtp.bcal = bcal; |
836 | if(triggerTypes[jj][6].size() > 0) trigger_tmp.gtp.en_thr = stoi(triggerTypes[jj][6]); |
837 | if((triggerTypes[jj].size() >= 12) && (triggerTypes[jj][11].size() > 0)) trigger_tmp.gtp.bcal_min = stoi(triggerTypes[jj][11]); |
838 | if((triggerTypes[jj].size() >= 13) && (triggerTypes[jj][12].size() > 0)) trigger_tmp.gtp.bcal_max = stoi(triggerTypes[jj][12]); |
839 | |
840 | |
841 | trig_found++; |
842 | } else { |
843 | cout << " Incorrect parameters for BFCAL trigger " << endl; |
844 | } |
845 | |
846 | } else if(triggerTypes[jj][1] == "ST"){ // ST |
847 | trigger_tmp.type = 0x4; |
848 | if(triggerTypes[jj][7].size() > 0) trigger_tmp.gtp.st_nhit = stoi(triggerTypes[jj][7]); |
849 | if((triggerTypes[jj].size() >= 14) && (triggerTypes[jj][13].size() > 0)) trigger_tmp.gtp.st_pattern = stoi(triggerTypes[jj][13]); |
850 | |
851 | trig_found++; |
852 | } else if(triggerTypes[jj][1] == "PS"){ // PS |
853 | trigger_tmp.type = 0x8; |
854 | if(triggerTypes[jj][7].size() > 0) trigger_tmp.gtp.ps_nhit = stoi(triggerTypes[jj][7]); |
855 | if((triggerTypes[jj].size() >= 14) && (triggerTypes[jj][13].size() > 0)) trigger_tmp.gtp.ps_pattern = stoi(triggerTypes[jj][13]); |
856 | trig_found++; |
857 | } else if(triggerTypes[jj][1] == "TAGH"){ // TAGH |
858 | trigger_tmp.type = 0x10; |
859 | if(triggerTypes[jj][7].size() > 0) trigger_tmp.gtp.tof_nhit = stoi(triggerTypes[jj][7]); |
860 | if((triggerTypes[jj].size() >= 14) && (triggerTypes[jj][13].size() > 0)){ |
861 | trigger_tmp.gtp.tagh_pattern = stoi(triggerTypes[jj][13],nullptr,0); |
862 | } |
863 | |
864 | trig_found++; |
865 | } else if(triggerTypes[jj][1] == "TOF"){ // TOF |
866 | trigger_tmp.type = 0x20; |
867 | if(triggerTypes[jj][7].size() > 0) trigger_tmp.gtp.ps_nhit = stoi(triggerTypes[jj][7]); |
868 | if((triggerTypes[jj].size() >= 14) && (triggerTypes[jj][13].size() > 0)) trigger_tmp.gtp.tof_pattern = stoi(triggerTypes[jj][13]); |
869 | trig_found++; |
870 | } else { |
871 | cout << " Incorrect Trigger type " << triggerTypes[jj][1] << endl; |
872 | } |
873 | |
874 | if(trig_found > 0){ |
875 | triggers_enabled.push_back(trigger_tmp); |
876 | } |
877 | |
878 | } // Trigger lane is enabled |
879 | } // Trigger types |
880 | |
881 | } |
882 | |
883 | |
884 | |
885 | // Load FCAL Trigger Masks |
886 | |
887 | string comDir = result.Sections["FCAL"].NameValues["FADC250_COM_DIR"]; |
888 | string comVer = result.Sections["FCAL"].NameValues["FADC250_COM_VER"]; |
889 | string userDir = result.Sections["FCAL"].NameValues["FADC250_USER_DIR"]; |
890 | string userVer = result.Sections["FCAL"].NameValues["FADC250_USER_VER"]; |
891 | |
892 | |
893 | for(int crate = 1; crate <= 12; crate++){ |
894 | |
895 | std::string s = std::to_string(crate); |
896 | |
897 | string comFileName = comDir + "/rocfcal" +s + "_fadc250_" + comVer + ".cnf"; |
898 | string userFileName = userDir + "/rocfcal" +s + "_" + userVer + ".cnf"; |
899 | |
900 | auto userFile = connection.GetFile(runnumber, userFileName); |
901 | |
902 | if(!userFile) { |
903 | // cout<<" USER File with name: "<< userFileName |
904 | // <<" doesn't exist (not associated) with run: "<< runnumber <<endl; |
905 | continue; |
906 | } |
907 | |
908 | |
909 | auto userParseResult = rcdb::ConfigParser::ParseWithSlots(userFile->GetContent(), "FADC250_SLOTS"); |
910 | |
911 | |
912 | for(unsigned int slot = 3; slot <= 21; slot++){ |
913 | |
914 | auto userValues = userParseResult.SectionsBySlotNumber[slot].NameVectors["FADC250_TRG_MASK"]; // Parse it and return |
915 | |
916 | if(userValues.size() > 0){ |
917 | for (unsigned int ch = 0; ch < userValues.size(); ++ch) { |
918 | |
919 | if(userValues[ch].size() == 0) continue; |
920 | |
921 | if(stoi(userValues[ch]) > 0){ |
922 | |
923 | uint32_t roc_id = 10 + crate; |
924 | |
925 | DTranslationTable::csc_t daq_index = {roc_id, slot, ch }; |
926 | |
927 | DTranslationTable::DChannelInfo channel_info; |
928 | |
929 | try { |
930 | channel_info = ttab[0]->GetDetectorIndex(daq_index); |
931 | } |
932 | |
933 | catch(...){ |
934 | if(print_messages) cout << "Exception: FCAL channel is not in the translation table " << " Crate = " << 10 + crate << " Slot = " << slot << |
935 | " Channel = " << ch << endl; |
936 | continue; |
937 | } |
938 | |
939 | fcal_mod tmp; |
940 | |
941 | tmp.roc = crate; |
942 | tmp.slot = slot; |
943 | tmp.ch = ch; |
944 | |
945 | tmp.row = channel_info.fcal.row; |
946 | tmp.col = channel_info.fcal.col; |
947 | |
948 | fcal_trig_mask.push_back(tmp); |
949 | } |
950 | |
951 | } |
952 | |
953 | } |
954 | |
955 | } // Loop over slots |
956 | } // Loop over crates |
957 | |
958 | |
959 | |
960 | |
961 | // Load BCAL Trigger Masks |
962 | |
963 | comDir = result.Sections["BCAL"].NameValues["FADC250_COM_DIR"]; |
964 | comVer = result.Sections["BCAL"].NameValues["FADC250_COM_VER"]; |
965 | userDir = result.Sections["BCAL"].NameValues["FADC250_USER_DIR"]; |
966 | userVer = result.Sections["BCAL"].NameValues["FADC250_USER_VER"]; |
967 | |
968 | |
969 | for(int crate = 1; crate <= 12; crate++){ |
970 | |
971 | if( (crate == 3) || (crate == 6) || (crate == 9) || (crate == 12)) continue; |
972 | |
973 | std::string s = std::to_string(crate); |
974 | |
975 | string comFileName = comDir + "/rocbcal" +s + "_fadc250_" + comVer + ".cnf"; |
976 | string userFileName = userDir + "/rocbcal" +s + "_" + userVer + ".cnf"; |
977 | |
978 | auto userFile = connection.GetFile(runnumber, userFileName); |
979 | |
980 | if(!userFile){ |
981 | // cout<<" USER File with name: "<< userFileName |
982 | // <<" doesn't exist (not associated) with run: "<< runnumber <<endl; |
983 | continue; |
984 | } |
985 | |
986 | |
987 | auto userParseResult = rcdb::ConfigParser::ParseWithSlots(userFile->GetContent(), "FADC250_SLOTS"); |
988 | |
989 | for(unsigned int slot = 3; slot <= 21; slot++){ |
990 | |
991 | auto userValues = userParseResult.SectionsBySlotNumber[slot].NameVectors["FADC250_TRG_MASK"]; // Parse it and return |
992 | |
993 | if(userValues.size() > 0){ |
994 | |
995 | for (unsigned int ch = 0; ch < userValues.size(); ++ch) { |
996 | |
997 | if(userValues[ch].size() == 0) continue; |
998 | |
999 | if(stoi(userValues[ch]) > 0){ |
1000 | |
1001 | uint32_t roc_id = 30 + crate; |
1002 | |
1003 | DTranslationTable::csc_t daq_index = {roc_id, slot, ch }; |
1004 | |
1005 | DTranslationTable::DChannelInfo channel_info; |
1006 | |
1007 | try { |
1008 | channel_info = ttab[0]->GetDetectorIndex(daq_index); |
1009 | } |
1010 | |
1011 | catch(...){ |
1012 | cout << "Exception: BCAL channel is not in the translation table " << " Crate = " << 30 + crate << " Slot = " << slot << |
1013 | " Channel = " << ch << endl; |
1014 | continue; |
1015 | } |
1016 | |
1017 | |
1018 | bcal_mod tmp; |
1019 | |
1020 | tmp.roc = crate; |
1021 | tmp.slot = slot; |
1022 | tmp.ch = ch; |
1023 | |
1024 | tmp.module = channel_info.bcal.module; |
1025 | tmp.layer = channel_info.bcal.layer; |
1026 | tmp.sector = channel_info.bcal.sector; |
1027 | tmp.end = channel_info.bcal.end; |
1028 | |
1029 | bcal_trig_mask.push_back(tmp); |
1030 | } |
1031 | |
1032 | } |
1033 | |
1034 | } |
1035 | |
1036 | } // Loop over slots |
1037 | } // Loop over crates |
1038 | |
1039 | return 0; |
1040 | |
1041 | #else // HAVE_RCDB |
1042 | |
1043 | return 10; // RCDB is not available |
1044 | |
1045 | #endif |
1046 | |
1047 | } |
1048 | |
1049 | |
1050 | int DL1MCTrigger_factory::SignalPulse(double en, double time, double amp_array[sample], int type){ |
1051 | |
1052 | |
1053 | // Parameterize and digitize pulse shapes. Sum up amplitudes |
1054 | // type = 1 - FCAL |
1055 | // type = 2 - BCAL |
1056 | |
1057 | float exp_par = 0.358; |
1058 | |
1059 | int pulse_length = 20; |
1060 | |
1061 | if(type == 2) exp_par = 0.18; |
1062 | |
1063 | int sample_first = (int)floor(time/time_stamp); |
1064 | |
1065 | // digitization range |
1066 | int ind_min = sample_first + 1; |
1067 | int ind_max = ind_min + pulse_length + 1; |
1068 | |
1069 | if( (ind_min > sample) || (ind_min < 0)){ |
1070 | // cout << " SignalPulse() FATAL error: time out of range " << time << " " << ind_min << endl; |
1071 | return 1; |
1072 | } |
1073 | |
1074 | if(ind_max > sample){ |
1075 | // cout << " SignalPulse(): ind_max set to maximum" << time << " " << ind_max << endl; |
1076 | ind_max = sample - 1; |
1077 | } |
1078 | |
1079 | for(int i = ind_min; i < ind_max; i++ ){ |
1080 | double adc_t = time_stamp*i - time; |
1081 | double amp = exp_par*exp_par*exp(-adc_t*exp_par)*adc_t; |
1082 | |
1083 | // amp_array[i] += (int)(amp*time_stamp*en + 0.5); |
1084 | // if(amp_array[i] > max_adc_bins){ |
1085 | // amp_array[i] = max_adc_bins; |
1086 | // } |
1087 | |
1088 | amp_array[i] += amp*time_stamp*en; |
1089 | |
1090 | } |
1091 | |
1092 | return 0; |
1093 | } |
1094 | |
1095 | int DL1MCTrigger_factory::GTP(int detector){ |
1096 | |
1097 | // 1 - FCAL |
1098 | // 2 - BCAL |
1099 | |
1100 | int INT_WINDOW = 20; |
1101 | |
1102 | switch(detector){ |
1103 | case 1: |
1104 | INT_WINDOW = FCAL_WINDOW; |
1105 | break; |
1106 | case 2: |
1107 | INT_WINDOW = BCAL_WINDOW; |
1108 | break; |
1109 | default: |
1110 | break; |
1111 | } |
1112 | |
1113 | int index_min = 0; |
1114 | int index_max = 0; |
1115 | |
1116 | for(unsigned int samp = 0; samp < sample; samp++){ |
1117 | |
1118 | index_max = samp; |
1119 | index_min = samp - INT_WINDOW; |
1120 | |
1121 | if(index_min < 0) index_min = 0; |
1122 | |
1123 | int energy_sum = 0; |
1124 | |
1125 | for(int ii = index_min; ii <= index_max; ii++){ |
1126 | if(detector == 1) |
1127 | energy_sum += fcal_ssp[ii]; |
1128 | else |
1129 | energy_sum += bcal_ssp[ii]; |
1130 | } |
1131 | |
1132 | if(detector == 1) |
1133 | fcal_gtp[samp] = energy_sum; |
1134 | else |
1135 | bcal_gtp[samp] = energy_sum; |
1136 | } |
1137 | |
1138 | return 0; |
1139 | |
1140 | } |
1141 | |
1142 | |
1143 | template <typename T> int DL1MCTrigger_factory::FADC_SSP(vector<T> merged_hits, int detector){ |
1144 | |
1145 | // 1 - FCAL |
1146 | // 2 - BCAL |
1147 | |
1148 | int EN_THR = 4096; |
1149 | int NSA = 10; |
1150 | int NSB = 3;; |
1151 | |
1152 | switch(detector){ |
1153 | case 1: |
1154 | EN_THR = FCAL_CELL_THR; |
1155 | NSA = FCAL_NSA; |
1156 | NSB = FCAL_NSB; |
1157 | break; |
1158 | case 2: |
1159 | EN_THR = BCAL_CELL_THR; |
1160 | NSA = BCAL_NSA; |
1161 | NSB = BCAL_NSB; |
1162 | break; |
1163 | default: |
1164 | break; |
1165 | } |
1166 | |
1167 | for(unsigned int hit = 0; hit < merged_hits.size(); hit++){ |
1168 | |
1169 | int index_min = -10; |
1170 | int index_max = -10; |
1171 | |
1172 | for(int ii = 0; ii < sample; ii++){ |
1173 | |
1174 | int pulse_found = 0; |
1175 | |
1176 | if(merged_hits[hit].adc_amp[ii] >= EN_THR){ |
1177 | pulse_found = 1; |
1178 | } else { |
1179 | continue; |
1180 | } |
1181 | |
1182 | index_min = ii - NSB; |
1183 | |
1184 | if(index_max > index_min) index_min = index_max + 1; |
1185 | |
1186 | index_max = ii + NSA - 1; |
1187 | |
1188 | if(index_min < 0) index_min = 0; |
1189 | |
1190 | if(index_max >= sample){ |
1191 | index_max = sample - 1; |
1192 | } |
1193 | |
1194 | int extend_nsa = 1; |
1195 | |
1196 | // Extend FADC range if needed |
1197 | |
1198 | while(extend_nsa){ |
1199 | |
1200 | int index_tmp = index_max + 1; |
1201 | |
1202 | if(index_tmp < sample){ |
1203 | if(merged_hits[hit].adc_amp[index_tmp] >= EN_THR){ |
1204 | index_max += NSA; |
1205 | } else extend_nsa = 0; |
1206 | } else extend_nsa = 0; |
1207 | } |
1208 | |
1209 | if(index_max >= sample) |
1210 | index_max = sample - 1; |
1211 | |
1212 | for(int kk = index_min; kk <= index_max; kk++){ |
1213 | if(detector == 1){ |
1214 | if((merged_hits[hit].adc_amp[kk] - 100) > 0) |
1215 | fcal_ssp[kk] += (merged_hits[hit].adc_amp[kk] - TRIG_BASELINE); |
1216 | } |
1217 | else if(detector == 2){ |
1218 | if((merged_hits[hit].adc_amp[kk] - 100) > 0) |
1219 | bcal_ssp[kk] += (merged_hits[hit].adc_amp[kk] - TRIG_BASELINE); |
1220 | } |
1221 | } |
1222 | |
1223 | if(pulse_found == 1){ |
1224 | ii = index_max + 1; |
1225 | pulse_found = 0; |
1226 | } |
1227 | |
1228 | } |
1229 | |
1230 | } |
1231 | |
1232 | return 0; |
1233 | } |
1234 | |
1235 | void DL1MCTrigger_factory::PrintTriggers(){ |
1236 | |
1237 | string detector; |
1238 | int nhit = 0; |
1239 | unsigned int pattern = 0; |
1240 | |
1241 | cout << endl << endl; |
1242 | cout << " ------------ Trigger Settings --------------- " << endl; |
1243 | cout << endl << endl; |
1244 | |
1245 | cout << "----------- FCAL ----------- " << endl << endl; |
1246 | |
1247 | cout << "FCAL_CELL_THR = " << setw(10) << FCAL_CELL_THR << endl; |
1248 | cout << "FCAL_NSA = " << setw(10) << FCAL_NSA << endl; |
1249 | cout << "FCAL_NSB = " << setw(10) << FCAL_NSB << endl; |
1250 | cout << "FCAL_WINDOW = " << setw(10) << FCAL_WINDOW << endl; |
1251 | |
1252 | cout << endl; |
1253 | |
1254 | cout << "----------- BCAL ----------- " << endl << endl; |
1255 | |
1256 | cout << "BCAL_CELL_THR = " << setw(10) << BCAL_CELL_THR << endl; |
1257 | cout << "BCAL_NSA = " << setw(10) << BCAL_NSA << endl; |
1258 | cout << "BCAL_NSB = " << setw(10) << BCAL_NSB << endl; |
1259 | cout << "BCAL_WINDOW = " << setw(10) << BCAL_WINDOW << endl; |
1260 | |
1261 | cout << endl << endl; |
1262 | |
1263 | |
1264 | if(triggers_enabled.size() > 0){ |
1265 | cout << "TYPE " << "FCAL_E " << "BCAL_E " << |
1266 | "EN_THR " << "NHIT " << "LANE " << "FCAL_EMIN " << "FCAL_EMAX " << |
1267 | "BCAL_EMIN " << "BCAL_EMAX " << "PATTERN " << endl; |
1268 | } |
1269 | |
1270 | |
1271 | for(unsigned int ii = 0; ii < triggers_enabled.size(); ii++){ |
1272 | |
1273 | switch(triggers_enabled[ii].type){ |
1274 | case 1: detector = "BFCAL "; |
1275 | break; |
1276 | case 2: detector = "BFCAL "; |
1277 | break; |
1278 | case 3: detector = "BFCAL "; |
1279 | break; |
1280 | case 4: detector = "ST "; |
1281 | nhit = triggers_enabled[ii].gtp.st_nhit; |
1282 | pattern = triggers_enabled[ii].gtp.st_pattern; |
1283 | break; |
1284 | case 8: detector = "PS "; |
1285 | nhit = triggers_enabled[ii].gtp.ps_nhit; |
1286 | pattern = triggers_enabled[ii].gtp.ps_pattern; |
1287 | break; |
1288 | case 16: detector = "TAGH "; |
1289 | nhit = triggers_enabled[ii].gtp.tagh_nhit; |
1290 | pattern = triggers_enabled[ii].gtp.tagh_pattern; |
1291 | break; |
1292 | case 32: detector = "TOF "; |
1293 | nhit = triggers_enabled[ii].gtp.tof_nhit; |
1294 | pattern = triggers_enabled[ii].gtp.tof_pattern; |
1295 | break; |
1296 | default: detector = "NONE "; |
1297 | cout << " Unknown detector ===== " << triggers_enabled[ii].type << endl; |
1298 | break; |
1299 | } |
1300 | |
1301 | cout << detector << setw(6) << |
1302 | triggers_enabled[ii].gtp.fcal << setw(9) << |
1303 | triggers_enabled[ii].gtp.bcal << setw(11) << |
1304 | triggers_enabled[ii].gtp.en_thr << setw(6) << |
1305 | nhit << setw(9) << |
1306 | triggers_enabled[ii].bit << setw(12) << |
1307 | triggers_enabled[ii].gtp.fcal_min << setw(14) << |
1308 | triggers_enabled[ii].gtp.fcal_max << setw(10) << |
1309 | triggers_enabled[ii].gtp.bcal_min << setw(14) << |
1310 | triggers_enabled[ii].gtp.bcal_max << setw(8) << |
1311 | hex << uppercase << "0x" << pattern << nouppercase << dec << endl; |
1312 | |
1313 | } |
1314 | |
1315 | cout << endl << endl; |
1316 | |
1317 | } |
1318 | |
1319 | |
1320 | |
1321 | int DL1MCTrigger_factory::FindTriggers(DL1MCTrigger *trigger){ |
1322 | |
1323 | int trig_found = 0; |
1324 | |
1325 | // Main production trigger |
1326 | for(unsigned int ii = 0; ii < triggers_enabled.size(); ii++){ |
1327 | |
1328 | if(triggers_enabled[ii].bit == 0){ // Main production trigger found |
1329 | |
1330 | int gtp_energy = 0; |
1331 | int bcal_energy = 0; |
1332 | |
1333 | for(unsigned int samp = 0; samp < sample; samp++){ |
1334 | |
1335 | int bcal_samp = samp - BCAL_OFFSET; |
1336 | |
1337 | if(bcal_samp < 0){ |
1338 | bcal_energy = 0; |
1339 | } else if(bcal_samp >= sample){ |
1340 | bcal_energy = 0; |
1341 | } else{ |
1342 | bcal_energy = bcal_gtp[bcal_samp]; |
1343 | } |
1344 | |
1345 | |
1346 | gtp_energy = triggers_enabled[ii].gtp.fcal*fcal_gtp[samp] + |
1347 | triggers_enabled[ii].gtp.bcal*bcal_energy; |
1348 | |
1349 | if(gtp_energy >= triggers_enabled[ii].gtp.en_thr){ |
1350 | |
1351 | if(triggers_enabled[ii].gtp.fcal_min > 0) { // FCAL > 0 |
1352 | if(fcal_gtp[samp] > triggers_enabled[ii].gtp.fcal_min){ |
1353 | trigger->trig_mask = (trigger->trig_mask | 0x1); |
1354 | trigger->trig_time[0] = samp - 25; |
1355 | trig_found++; |
1356 | } |
1357 | } else { |
1358 | |
1359 | trigger->trig_mask = (trigger->trig_mask | 0x1); |
1360 | trigger->trig_time[0] = samp - 25; |
1361 | trig_found++; |
1362 | } |
1363 | |
1364 | break; |
1365 | |
1366 | } // Check energy threshold |
1367 | } |
1368 | } // Trigger Bit 0 |
1369 | |
1370 | |
1371 | if(triggers_enabled[ii].bit == 2){ // BCAL trigger found |
1372 | |
1373 | int gtp_energy = 0; |
1374 | int bcal_energy = 0; |
1375 | |
1376 | for(unsigned int samp = 0; samp < sample; samp++){ |
1377 | |
1378 | int bcal_samp = samp - BCAL_OFFSET; |
1379 | |
1380 | if(bcal_samp < 0){ |
1381 | bcal_energy = 0; |
1382 | } else if(bcal_samp >= sample){ |
1383 | bcal_energy = 0; |
1384 | } else{ |
1385 | bcal_energy = bcal_gtp[bcal_samp]; |
1386 | } |
1387 | |
1388 | gtp_energy = triggers_enabled[ii].gtp.bcal*bcal_energy; |
1389 | |
1390 | if(gtp_energy >= triggers_enabled[ii].gtp.en_thr){ |
1391 | |
1392 | trigger->trig_mask = (trigger->trig_mask | 0x4); |
1393 | trigger->trig_time[2] = samp - 25; |
1394 | trig_found++; |
1395 | |
1396 | break; |
1397 | |
1398 | } // Energy threshold |
1399 | } |
1400 | } // Trigger Bit 3 |
1401 | |
1402 | } // Loop over triggers found in the config file |
1403 | |
1404 | return trig_found; |
1405 | |
1406 | } |
1407 | |
1408 | |
1409 | // Fill fcal calibration tables similar to FCALHit factory |
1410 | void DL1MCTrigger_factory::LoadFCALConst(fcal_constants_t &table, const vector<double> &fcal_const_ch, |
1411 | const DFCALGeometry &fcalGeom){ |
1412 | for (int ch = 0; ch < static_cast<int>(fcal_const_ch.size()); ch++) { |
1413 | int row = fcalGeom.row(ch); |
1414 | int col = fcalGeom.column(ch); |
1415 | table[row][col] = fcal_const_ch[ch]; |
1416 | } |
1417 | |
1418 | } |
1419 | |
1420 | void DL1MCTrigger_factory::Digitize(double adc_amp[sample], int adc_count[sample]){ |
1421 | |
1422 | for(int samp = 0; samp < sample; samp++ ){ |
1423 | |
1424 | adc_count[samp] += (int)(adc_amp[samp] + TRIG_BASELINE + 0.5); |
1425 | |
1426 | if(adc_count[samp] > max_adc_bins) |
1427 | adc_count[samp] = max_adc_bins; |
1428 | |
1429 | } |
1430 | } |
1431 | |
1432 | |
1433 | void DL1MCTrigger_factory::AddBaseline(double adc_amp[sample], double pedestal, DRandom2 &gDRandom){ |
1434 | |
1435 | double pedestal_correct = pedestal - TRIG_BASELINE; |
1436 | |
1437 | for(int samp = 0; samp < sample; samp++ ){ |
1438 | double tmp = gDRandom.Gaus(pedestal_correct, pedestal_sigma); |
1439 | adc_amp[samp] += tmp; |
1440 | // cout << " " << tmp; |
1441 | } |
1442 | |
1443 | if(debug){ |
1444 | cout << endl; |
1445 | cout << " Corrected pedestals = " << pedestal_correct << " " << adc_amp[sample - 2] |
1446 | << " " << pedestal_sigma << endl; |
1447 | } |
1448 | |
1449 | } |
1450 | |
1451 | |
1452 | |
1453 | void DL1MCTrigger_factory::GetSeeds(JEventLoop *loop, uint64_t eventnumber, UInt_t &seed1, UInt_t &seed2, UInt_t &seed3){ |
1454 | |
1455 | // Use seeds similar to mcsmear |
1456 | |
1457 | JEvent& event = loop->GetJEvent(); |
1458 | |
1459 | JEventSource *source = event.GetJEventSource(); |
1460 | |
1461 | DEventSourceHDDM *hddm_source = dynamic_cast<DEventSourceHDDM*>(source); |
1462 | |
1463 | if (!hddm_source) { |
1464 | |
1465 | cerr << "DL1MCTrigger_factory: This program MUST be used with an HDDM file as input!" << endl; |
1466 | cerr << " Default seeds will be used for the random generator " << endl; |
1467 | seed1 = 259921049 + eventnumber; |
1468 | seed2 = 442249570 + eventnumber; |
1469 | seed3 = 709975946 + eventnumber; |
1470 | } else { |
1471 | |
1472 | hddm_s::HDDM *record = (hddm_s::HDDM*)event.GetRef(); |
1473 | if (!record){ |
1474 | seed1 = 259921049 + eventnumber; |
1475 | seed2 = 442249570 + eventnumber; |
1476 | seed3 = 709975946 + eventnumber; |
1477 | } else { |
1478 | |
1479 | |
1480 | hddm_s::ReactionList::iterator reiter = record->getReactions().begin(); |
1481 | |
1482 | hddm_s::Random my_rand = reiter->getRandom(); |
1483 | |
1484 | // Copy seeds from event record to local variables |
1485 | seed1 = my_rand.getSeed1(); |
1486 | seed2 = my_rand.getSeed2(); |
1487 | seed3 = my_rand.getSeed3(); |
1488 | |
1489 | // If no seeds are stored in the hddm file, generate them in the same way |
1490 | // as in mcsmear |
1491 | if ((seed1 == 0) || (seed2 == 0) || (seed3 == 0)){ |
1492 | uint64_t eventNo = record->getPhysicsEvent().getEventNo(); |
1493 | seed1 = 259921049 + eventNo; |
1494 | seed2 = 442249570 + eventNo; |
1495 | seed3 = 709975946 + eventNo; |
1496 | |
1497 | } |
1498 | |
1499 | } // Record doesn't exist |
1500 | } // Not an HDDM file |
1501 | } |