Bug Summary

File:libraries/TRIGGER/DL1MCTrigger_factory.cc
Location:line 560, column 6
Description:Value stored to 'status' is never read

Annotated Source Code

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