Bug Summary

File:/volatile/halld/gluex/nightly/2024-08-19/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/TRIGGER/DL1MCTrigger_factory_DATA.cc
Location:line 637, column 7
Description:Value stored to 'status' is never read

Annotated Source Code

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