Bug Summary

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