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