« back to all changes in this revision
Viewing changes to src/map/MapCppCkovMCDigitizer/MapCppCkovMCDigitizer.cc
- browse files at revision 697.93.1
- compare with another revision
- download diff
- download tarball
- view history from revision 697.93.1
- Committer: Ao Liu
- Date: 2015-09-01 03:14:59 UTC
- mfrom: (697.90.8 trunk)
- mto: (697.90.11 trunk)
- mto: This revision was merged to the branch mainline in revision 711.
- Revision ID: aoliu@fnal.gov-20150901031459-i0stzz3bzo7b06zh
CkovMCDigitizer with the MAUS::Data and better handling of the hits
- files added:
- src/map/MapCppCkovMCDigitizer/ckov_geom.dat
- src/map/MapCppMCReconSetup
- files removed:
- src/map/MapCppTOFMCDigitizer/mc_test.dat
- third_party/install/lib64
- third_party/install/lib64/Geant4-9.6.2
- third_party/install/lib64/Geant4-9.6.2/Linux-g++
- third_party/install/lib64/Geant4-9.6.2/Modules
- tmp000
- files modified:
- .bzrignore
added
removed
src/map/MapCppCkovMCDigitizer/MapCppCkovMCDigitizer.cc
22
22
#include "API/PyWrapMapBase.hh"
23
23
#include "src/common_cpp/DataStructure/RunHeaderData.hh"
24
24
#include "src/common_cpp/DataStructure/RunHeader.hh"
25
#include <numeric>
25
26
26
27
namespace MAUS {
27
28
PyMODINIT_FUNC init_MapCppCkovMCDigitizer(void) {
30
31
}
31
32
32
33
MapCppCkovMCDigitizer::MapCppCkovMCDigitizer()
33
: MapBase<Json::Value>("MapCppCkovMCDigitizer") {
34
: MapBase<Data>("MapCppCkovMCDigitizer") {
34
35
}
35
36
36
37
//////////////////////////////////////////////////////////////////////
39
40
// Before we go on we generate a new TRandom3;
40
41
rnd = new TRandom3;
41
42
// print level
42
fDebug = false;
43
fDebug = true;
43
44
44
45
Json::Reader reader;
45
46
66
67
geo_module = new MiceModule(filename);
67
68
ckov_modules = geo_module->findModulesByPropertyString("SensitiveDetector", "CKOV");
68
69
69
_station_index.push_back("A");
70
_station_index.push_back("B");
70
// initialize thesholds and conversion factors
71
// Ckov A
72
muon_thrhold[0] = 275.0;
73
charge_per_pe[0] = 17.3;
74
// Ckov B
75
muon_thrhold[1] = 210.0;
76
charge_per_pe[1] = 26;
77
78
// scaling factor to data
79
scaling = 0.935;
80
81
// npe->adc conversion factor
82
ckovNpeToAdc = 23.0;
71
83
}
72
84
73
85
//////////////////////////////////////////////////////////////////////
75
87
}
76
88
77
89
//////////////////////////////////////////////////////////////////////
78
void MapCppCkovMCDigitizer::_process(Json::Value* document) const {
90
void MapCppCkovMCDigitizer::_process(MAUS::Data* data) const {
79
91
80
Json::Value& root = *document;
81
// check sanity of json input file and mc brach
82
Json::Value mc = check_sanity_mc(*document);
83
92
84
93
// --- Below are from TOF:
85
94
// all_tof_digits store all the tof hits
87
96
// ---
88
97
// For the Ckov, it's similar. Just change TOF to Ckov.
89
98
//
90
std::vector<Json::Value> all_ckov_digits;
91
all_ckov_digits.push_back(Json::Value());
99
100
// Get spill, break if there's no DAQ data
101
Spill *spill = data->GetSpill();
102
if (spill->GetMCEvents() == NULL)
103
return;
104
105
MCEventPArray *mcEvts = spill->GetMCEvents();
106
int nPartEvents = spill->GetMCEventSize();
107
if (nPartEvents == 0)
108
return;
109
110
ReconEventPArray *recEvts = spill->GetReconEvents();
111
112
// Resize the recon event to hold mc events
113
if (spill->GetReconEventSize() == 0) { // No recEvts yet
114
for (int iPe = 0; iPe < nPartEvents; iPe++) {
115
recEvts->push_back(new ReconEvent);
116
}
117
}
118
119
if (fDebug) std::cerr << "mc numevt = i" << mcEvts->size() << std::endl;
92
120
93
121
// loop over events
94
if (fDebug) printf("mc numevt = %i", mc.size());
95
for ( unsigned int i = 0; i < mc.size(); i++ ) {
96
Json::Value particle = mc[i];
97
double gentime = particle["primary"]["time"].asDouble();
98
if (fDebug) std::cout << "evt: " << i << " pcle= " << particle << std::endl;
99
100
// hits for this event
101
Json::Value _hits = particle["ckov_hits"];
102
// if _hits.size() = 0, the make_digits and fill_ckov_evt functions will
103
// return null. The Ckov recon expects something - either null or data, can't just skip
104
all_ckov_digits.clear();
105
106
// pick out tof hits, digitize and store
107
all_ckov_digits = make_ckov_digits(_hits, gentime, *document);
108
109
// Does Ckov need this pixel thing? Maybe not.
110
// std::string strig = findTriggerPixel(all_tof_digits);
111
112
// ------>NOW, go through the stations and fill Ckov events
113
// real data tof digits look like so:
114
// "digits":{"tof1":[ [evt0..{pmt0},{pmt1}..],[evt1..]],"tof0":[[evt0]..],tof2..}
115
// loop over ckov stations
116
Json::Value ckov_digitAB(Json::arrayValue);
117
ckov_digitAB.clear();
118
for (int snum = 0; snum < 2; ++snum) {
119
for (unsigned int kk = 0; kk < all_ckov_digits.size(); ++kk) {
120
// check that this digit belongs to the station we are trying to fill
121
if (all_ckov_digits[kk]["station"].asInt() != snum) continue;
122
ckov_digitAB.append(fill_ckov_evt(i, snum, all_ckov_digits, kk));
123
}
124
} // end loop over stations
125
root["recon_events"][i]["ckov_event"]["ckov_digits"] = ckov_digitAB;
122
for ( unsigned int i = 0; i < mcEvts->size(); i++ ) {
123
// Get Ckov hits for this event
124
CkovHitArray *hit_array = spill->GetAnMCEvent(i).GetCkovHits();
125
if (fDebug) printf("\nNumber of hits: %i\n", hit_array->size());
126
127
double gentime = spill->GetAnMCEvent(i).GetPrimary()->GetTime();
128
CkovEvent *evt = new CkovEvent();
129
(*recEvts)[i]->SetCkovEvent(evt);
130
131
// process only if there are hits
132
if (hit_array) {
133
134
// pick out ckov hits, digitize and store
135
multi_ckov_dig all_ckov_dig = make_ckov_digits(hit_array, gentime);
136
137
CkovDigitArray* digit_array = evt->GetCkovDigitArrayPtr();
138
fill_ckov_evt(i, all_ckov_dig, digit_array);
139
140
} // end check-if-hits
126
141
} // end loop over events
127
142
}
128
143
129
144
//////////////////////////////////////////////////////////////////////
130
std::vector<Json::Value>
131
MapCppCkovMCDigitizer::make_ckov_digits(Json::Value hits, double gentime,
132
Json::Value& root) const {
133
std::vector<Json::Value> ckov_digits;
134
ckov_digits.clear();
135
if (hits.size() == 0) return ckov_digits;
136
137
for (unsigned int j = 0; j < hits.size(); ++j) { // j-th hit
138
Json::Value hit = hits[j];
139
145
multi_ckov_dig MapCppCkovMCDigitizer::make_ckov_digits(CkovHitArray* hit_array,
146
double gentime) const {
147
148
multi_ckov_dig all_ckov_dig(0);
149
if (hit_array->size() == 0) return all_ckov_dig;
150
151
for (unsigned int j = 0; j < hit_array->size(); ++j) { // j-th hit
152
CkovHit hit = hit_array->at(j);
140
153
if (fDebug) std::cout << "=================== hit# " << j << std::endl;
141
154
142
// make sure we can get the station/slab info
143
if (!hit.isMember("channel_id"))
155
// make sure we can get the station number
156
if (!hit.GetChannelId())
144
157
throw(Exception(Exception::recoverable,
145
158
"No channel_id in hit",
146
159
"MapCppCkovMCDigitizer::make_ckov_digits"));
147
if (!hit.isMember("momentum"))
148
throw(Exception(Exception::recoverable,
149
"No momentum in hit",
150
"MapCppCkovMCDigitizer::make_ckov_digits"));
151
if (!hit.isMember("time"))
152
throw(Exception(Exception::recoverable,
153
"No time in hit",
154
"MapCppCkovMCDigitizer::make_ckov_digits"));
155
Json::Value channel_id = hit["channel_id"];
156
160
157
if (!hit.isMember("energy_deposited"))
158
throw(Exception(Exception::recoverable,
159
"No energy_deposited in hit",
160
"MapCppCkovMCDigitizer::make_ckov_digits"));
161
double edep = hit["energy_deposited"].asDouble();
161
double edep = hit.GetEnergyDeposited();
162
162
163
163
// Ckov hits in the CkovChannelIdProcessor of Json Parser registers a property
164
164
// "station_number" in channel_id
165
int stn = hit["channel_id"]["station_number"].asInt();
165
int stn = hit.GetChannelId()->GetStation();
166
166
167
167
// find the geo module corresponding to this hit
168
168
const MiceModule* hit_module = NULL;
179
179
"MapCppCkovMCDigitizer::make_ckov_digits"));
180
180
181
181
// now get the position of the hit
182
Hep3Vector hpos = JsonWrapper::JsonToThreeVector(hit["position"]);
182
ThreeVector tpos = hit.GetPosition();
183
Hep3Vector hpos(tpos.x(), tpos.y(), tpos.z());
183
184
184
185
// get the pmt positions from the geometry
185
186
// pmt positions
194
195
pmtName = "PMT2Pos";
195
196
Hep3Vector pmt2Pos = hit_module->propertyHep3Vector(pmtName);
196
197
// pmt 3
197
pmtName = "PMT2Pos";
198
pmtName = "PMT3Pos";
198
199
Hep3Vector pmt3Pos = hit_module->propertyHep3Vector(pmtName);
199
200
//
200
201
206
207
207
208
// Here we uniformly distribute all the pe to the 4 pmts;
208
209
// The distance is used to calculate the time of the signal.
209
printf("Passing hit number %i", j);
210
if (fDebug) printf("Passing hit number %i", j);
210
211
double npe = get_npe(hit);
211
212
double npe0, npe1, npe2, npe3;
212
213
npe0 = npe1 = npe2 = npe3 = npe/4;
213
214
214
215
// get the hit time, speed of c in the ckov, and define the pmt resolution;
215
double htime = hit["time"].asDouble() - gentime;
216
double htime = hit.GetTime() - gentime;
216
217
double csp = (stn == 0 ? 3.0e8/1.069 : 3.0e8/1.112);
217
218
double pmt_res = 0.1;
218
219
225
226
double tdc_factor = 0.025;
226
227
227
228
// convert to tdc
229
// I believe the ckov only has caen v1731 flash adc's
228
230
int tdc0 = static_cast<int>(time0 / tdc_factor);
229
231
int tdc1 = static_cast<int>(time1 / tdc_factor);
230
232
int tdc2 = static_cast<int>(time2 / tdc_factor);
234
236
<< " " << tdc2 << " " << tdc3 << std::endl;
235
237
}
236
238
237
Json::Value tmpDigit;
238
tmpDigit.clear();
239
tmpDigit["channel_id"] = hit["channel_id"];
240
tmpDigit["momentum"] = hit["momentum"];
241
tmpDigit["time"] = hit["time"];
242
tmpDigit["position"]=hit["position"];
243
tmpDigit["isUsed"] = 0;
244
tmpDigit["npe"] = npe;
245
246
// set the key for this channel
247
tmpDigit["station"] = stn;
248
tmpDigit["npe0"] = npe0;
249
tmpDigit["leading_time0"] = tdc0;
250
tmpDigit["raw_time0"] = time0;
251
tmpDigit["npe1"] = npe1;
252
tmpDigit["leading_time1"] = tdc1;
253
tmpDigit["raw_time1"] = time1;
254
tmpDigit["npe2"] = npe2;
255
tmpDigit["leading_time2"] = tdc2;
256
tmpDigit["raw_time2"] = time2;
257
tmpDigit["npe3"] = npe3;
258
tmpDigit["leading_time3"] = tdc3;
259
tmpDigit["raw_time3"] = time3;
260
ckov_digits.push_back(tmpDigit);
261
} // for (unsigned int j = 0; j < hits.size(); ++j)
262
263
return ckov_digits;
264
}
265
266
//////////////////////////////////////////////////////////////////////
267
Json::Value MapCppCkovMCDigitizer::check_sanity_mc(Json::Value& root) const {
268
// Check if the JSON document can be parsed, else return error only
269
// Check if the JSON document has a 'mc' branch, else return error
270
if (!root.isMember("mc_events")) {
271
throw MAUS::Exception(Exception::recoverable,
272
"Could not find an MC branch to simulate.",
273
"MapCppCkovMCDigitizer::check_sanity_mc");
274
}
275
276
Json::Value mc = root.get("mc_events", 0);
277
// Check if JSON document is of the right type, else return error
278
if (!mc.isArray()) {
279
throw MAUS::Exception(Exception::recoverable,
280
"MC branch was not an array.",
281
"MapCppCkovMCDigitizer::check_sanity_mc");
282
}
283
return mc;
284
}
285
286
//////////////////////////////////////////////////////////////////////
287
double MapCppCkovMCDigitizer::get_npe(Json::Value hit) const {
239
one_ckov_dig a_ckov_dig;
240
a_ckov_dig.fChannelId = hit.GetChannelId();
241
a_ckov_dig.fMom = hit.GetMomentum();
242
a_ckov_dig.fTime = hit.GetTime();
243
a_ckov_dig.fPos = hit.GetPosition();
244
a_ckov_dig.fIsUsed = 0;
245
a_ckov_dig.fNpe = npe;
246
247
// set the key for this channel
248
a_ckov_dig.fStation = stn;
249
a_ckov_dig.fNpe0 = npe0;
250
a_ckov_dig.fLeadingTime0 = tdc0;
251
a_ckov_dig.fRawTime0 = time0;
252
a_ckov_dig.fNpe1 = npe1;
253
a_ckov_dig.fLeadingTime1 = tdc1;
254
a_ckov_dig.fRawTime1 = time1;
255
a_ckov_dig.fNpe2 = npe2;
256
a_ckov_dig.fLeadingTime2 = tdc2;
257
a_ckov_dig.fRawTime2 = time2;
258
a_ckov_dig.fNpe3 = npe3;
259
a_ckov_dig.fLeadingTime3 = tdc3;
260
a_ckov_dig.fRawTime3 = time3;
261
all_ckov_dig.push_back(a_ckov_dig);
262
} // for (unsigned int j = 0; j < hit_array.size(); ++j)
263
264
return all_ckov_dig;
265
}
266
//////////////////////////////////////////////////////////////////////
267
double MapCppCkovMCDigitizer::get_npe(CkovHit& hit) const {
288
268
289
269
double nphot = 0.;
290
270
double nptmp = 0.;
291
double muon_thrhold = 0.;
292
double charge_per_pe = 0.;
293
// scaling factor to data
294
double scaling = 0.935;
295
271
// mass of the hit particle;
296
double particle_mass = hit["mass"].asDouble();
272
double particle_mass = hit.GetMass();
297
273
// momentum of the hit particle;
298
double px = hit["momentum"]["x"].asDouble();
299
double py = hit["momentum"]["y"].asDouble();
300
double pz = hit["momentum"]["z"].asDouble();
274
double px = hit.GetMomentum().x();
275
double py = hit.GetMomentum().y();
276
double pz = hit.GetMomentum().z();
301
277
double p_tot = sqrt(px*px + py*py + pz*pz);
302
278
303
if (hit["channel_id"]["station_number"].asInt() == 0) {
304
muon_thrhold = 275.0;
305
charge_per_pe = 17.3;
306
} else if (hit["channel_id"]["station_number"].asInt() == 1) {
307
muon_thrhold = 210.0;
308
charge_per_pe = 26;
309
}
279
int hitStation = hit.GetChannelId()->GetStation();
280
int hitPid = hit.GetParticleId();
310
281
311
282
// momentum threshold, see Lucien's Python code.
312
double pp_threshold = muon_thrhold * particle_mass / 105.6;
283
double pp_threshold = muon_thrhold[hitStation] * particle_mass / 105.6;
313
284
314
285
// smearing constant
315
286
const double smear_const = 0.5;
325
296
if (p_tot >= pp_threshold) {
326
297
327
298
// Above threshold;
328
double pred = scaling * charge_per_pe *
299
double pred = scaling * charge_per_pe[hitStation] *
329
300
(1.0 - (pp_threshold/p_tot)*(pp_threshold/p_tot)) + 0.5;
330
301
double pred_single = pred/4;
331
302
333
304
rms_sum = 0;
334
305
335
306
// Electron shower:
336
if (hit["particle_id"].asInt() == 11 || hit["particle_id"].asInt() == -11) {
307
if (hitPid == 11 || hitPid == -11) {
337
308
338
309
double rnd_number = rnd->Rndm();
339
if (hit["channel_id"]["station_number"].asInt() == 0) {
310
if (hitStation == 0) {
340
311
if (rnd_number <= 0.03)
341
312
pred *= 2;
342
313
else if (rnd_number <= 0.039 && rnd_number > 0.03)
365
336
rms = sqrt(rms_sum);
366
337
rms = (rms < 0.2 ? 0.2 : rms);
367
338
npssn = npssn_sum;
368
if (hit["particle_id"].asInt() == 13 || hit["particle_id"].asInt() == -13 ||
369
hit["particle_id"].asInt() == 211 || hit["particle_id"].asInt() == -211) {
339
if (hitPid == 13 || hitPid == -13 ||
340
hitPid == 211 || hitPid == -211) {
370
341
if (rnd->Rndm() < 0.06) {
371
342
npssn -= log(rnd->Rndm())/0.2;
372
343
}
404
375
}
405
376
406
377
//////////////////////////////////////////////////////////////////////
407
Json::Value MapCppCkovMCDigitizer::fill_ckov_evt(int evnum,
408
int snum,
409
std::vector<Json::Value> all_ckov_digits,
410
int i)
411
const {
412
Json::Value ckov_digits;
413
ckov_digits.clear();
414
415
// return null if this evt had no ckov hits
416
if (all_ckov_digits.size() == 0) return ckov_digits;
417
418
double npe;
419
int ndigs = 0;
420
421
Json::Value digitA;
422
Json::Value digitB;
423
digitA.clear();
424
digitB.clear();
425
426
// check if the digit has been used
427
if (all_ckov_digits[i]["isUsed"].asInt() == 0) {
428
429
ndigs++;
430
431
npe = all_ckov_digits[i]["npe"].asDouble();
432
// 23 is the ADC factor.
433
int adc = static_cast<int>(npe * 23);
434
435
// NOTE: needs tweaking/verifying -- DR 3/15
436
// NOTE: if tof needs tweaking, Ckov also needs it. -- frankliuao 8/15
437
// Initialize digits for both stations:
438
digitA["position_min_0"] = 0;
439
digitA["position_min_1"] = 0;
440
digitA["position_min_2"] = 0;
441
digitA["position_min_3"] = 0;
442
digitA["pulse_0"] = 0;
443
digitA["pulse_1"] = 0;
444
digitA["pulse_2"] = 0;
445
digitA["pulse_3"] = 0;
446
digitA["coincidences"] = 0;
447
digitA["total_charge"] = 0;
448
digitA["arrival_time_0"] = 0;
449
digitA["arrival_time_1"] = 0;
450
digitA["arrival_time_2"] = 0;
451
digitA["arrival_time_3"] = 0;
452
digitA["part_event_number"] = 0;
453
digitA["number_of_pes"] = 0;
454
455
456
digitB["position_min_4"] = 0;
457
digitB["position_min_5"] = 0;
458
digitB["position_min_6"] = 0;
459
digitB["position_min_7"] = 0;
460
digitB["pulse_4"] = 0;
461
digitB["pulse_5"] = 0;
462
digitB["pulse_6"] = 0;
463
digitB["pulse_7"] = 0;
464
digitB["total_charge"] = 0;
465
digitB["coincidences"] = 0;
466
digitB["arrival_time_4"] = 0;
467
digitB["arrival_time_5"] = 0;
468
digitB["arrival_time_6"] = 0;
469
digitB["arrival_time_7"] = 0;
470
digitB["part_event_number"] = 0;
471
digitB["number_of_pes"] = 0;
472
473
474
if (snum == 0) {
475
digitA["arrival_time_0"] = all_ckov_digits[i]["leading_time0"];
476
digitA["arrival_time_1"] = all_ckov_digits[i]["leading_time1"];
477
digitA["arrival_time_2"] = all_ckov_digits[i]["leading_time2"];
478
digitA["arrival_time_3"] = all_ckov_digits[i]["leading_time3"];
479
digitA["total_charge"] = adc;
480
digitA["part_event_number"] = evnum;
481
digitA["number_of_pes"] = npe;
482
} else {
483
digitB["arrival_time_4"] = all_ckov_digits[i]["leading_time0"];
484
digitB["arrival_time_5"] = all_ckov_digits[i]["leading_time1"];
485
digitB["arrival_time_6"] = all_ckov_digits[i]["leading_time2"];
486
digitB["arrival_time_7"] = all_ckov_digits[i]["leading_time3"];
487
digitB["total_charge"] = adc;
488
digitB["part_event_number"] = evnum;
489
digitB["number_of_pes"] = npe;
490
}
491
492
all_ckov_digits[i]["isUsed"] = 1;
493
} // #if (all_ckov_digits[i]["isUsed"].asInt() == 0)
494
ckov_digits["A"] = digitA;
495
ckov_digits["B"] = digitB;
496
return ckov_digits;
378
void MapCppCkovMCDigitizer::fill_ckov_evt(int evnum,
379
multi_ckov_dig& all_ckov_dig,
380
CkovDigitArray* digit_array) const {
381
// return null if this evt had no ckov hits, no need to go further;
382
if (all_ckov_dig.size() == 0) return;
383
384
/* If there are hits:
385
* For this given event number, there might be multiple hits;
386
* If the hits belong to 2 stations, they should be filled simultaneously to
387
* the digits, but in their respective structure;
388
* Otherwise, we need to integrate the digits together and get an average;
389
*/
390
std::vector<double> npe_each_hit_a;
391
std::vector<double> npe_each_hit_b;
392
std::vector<double> arrival_time0_each_hit_a;
393
std::vector<double> arrival_time1_each_hit_a;
394
std::vector<double> arrival_time2_each_hit_a;
395
std::vector<double> arrival_time3_each_hit_a;
396
std::vector<double> arrival_time4_each_hit_b;
397
std::vector<double> arrival_time5_each_hit_b;
398
std::vector<double> arrival_time6_each_hit_b;
399
std::vector<double> arrival_time7_each_hit_b;
400
401
// Count how many hits each station has
402
int num_of_hits_a = 0;
403
int num_of_hits_b = 0;
404
405
for (int snum = 0; snum < 2; ++snum) {
406
for (unsigned int kk = 0; kk < all_ckov_dig.size(); ++kk) {
407
// check that this hit has not already been used/filled
408
if (all_ckov_dig[kk].fStation != snum) continue;
409
if (all_ckov_dig[kk].fIsUsed == 0) {
410
// 23 is the ADC factor.
411
// the factor - ckovNpeToAdc is set at birth, defined in header
412
// keeping for-now-hardcoded definitions in one place
413
// -- DR 2015-08-31
414
415
if (snum == 0) {
416
arrival_time0_each_hit_a.push_back(all_ckov_dig[kk].fLeadingTime0);
417
arrival_time1_each_hit_a.push_back(all_ckov_dig[kk].fLeadingTime1);
418
arrival_time2_each_hit_a.push_back(all_ckov_dig[kk].fLeadingTime2);
419
arrival_time3_each_hit_a.push_back(all_ckov_dig[kk].fLeadingTime3);
420
npe_each_hit_a.push_back(all_ckov_dig[kk].fNpe);
421
num_of_hits_a++;
422
} else {
423
arrival_time4_each_hit_b.push_back(all_ckov_dig[kk].fLeadingTime0);
424
arrival_time5_each_hit_b.push_back(all_ckov_dig[kk].fLeadingTime1);
425
arrival_time6_each_hit_b.push_back(all_ckov_dig[kk].fLeadingTime2);
426
arrival_time7_each_hit_b.push_back(all_ckov_dig[kk].fLeadingTime3);
427
npe_each_hit_b.push_back(all_ckov_dig[kk].fNpe);
428
num_of_hits_b++;
429
}
430
all_ckov_dig[kk].fIsUsed = true;
431
} // end check if-digit-used
432
} // end loop over tmp digits
433
} // end loop over stations
434
435
// Now, if either num_of_hits_a(b) == 0, write the default values to the station:
436
// Can not be 0 at the same time. If there are no hits, it won't reach this far.
437
CkovA digitA;
438
CkovB digitB;
439
if (num_of_hits_a == 0) {
440
441
digitA.SetArrivalTime0(-999);
442
digitA.SetArrivalTime1(-999);
443
digitA.SetArrivalTime2(-999);
444
digitA.SetArrivalTime3(-999);
445
digitA.SetTotalCharge(-999);
446
digitA.SetPartEventNumber(evnum);
447
digitA.SetNumberOfPes(-999);
448
digitA.SetPositionMin0(-999);
449
digitA.SetPositionMin1(-999);
450
digitA.SetPositionMin2(-999);
451
digitA.SetPositionMin3(-999);
452
digitA.SetPulse0(-999);
453
digitA.SetPulse1(-999);
454
digitA.SetPulse2(-999);
455
digitA.SetPulse3(-999);
456
digitA.SetCoincidences(-999);
457
458
digitB.SetArrivalTime4(std::accumulate(arrival_time4_each_hit_b.begin(),
459
arrival_time4_each_hit_b.end(),
460
0.0) / arrival_time4_each_hit_b.size());
461
digitB.SetArrivalTime5(std::accumulate(arrival_time5_each_hit_b.begin(),
462
arrival_time5_each_hit_b.end(),
463
0.0) / arrival_time5_each_hit_b.size());
464
digitB.SetArrivalTime6(std::accumulate(arrival_time6_each_hit_b.begin(),
465
arrival_time6_each_hit_b.end(),
466
0.0) / arrival_time6_each_hit_b.size());
467
digitB.SetArrivalTime7(std::accumulate(arrival_time7_each_hit_b.begin(),
468
arrival_time7_each_hit_b.end(),
469
0.0) / arrival_time7_each_hit_b.size());
470
double total_npe = std::accumulate(npe_each_hit_b.begin(),
471
npe_each_hit_b.end(),
472
0.0);
473
digitB.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
474
digitB.SetPartEventNumber(evnum);
475
digitB.SetNumberOfPes(total_npe);
476
digitB.SetPositionMin4(0);
477
digitB.SetPositionMin5(0);
478
digitB.SetPositionMin6(0);
479
digitB.SetPositionMin7(0);
480
digitB.SetPulse4(0);
481
digitB.SetPulse5(0);
482
digitB.SetPulse6(0);
483
digitB.SetPulse7(0);
484
digitB.SetCoincidences(0);
485
} else if (num_of_hits_b == 0) {
486
digitB.SetArrivalTime4(-999);
487
digitB.SetArrivalTime5(-999);
488
digitB.SetArrivalTime6(-999);
489
digitB.SetArrivalTime7(-999);
490
digitB.SetTotalCharge(-999);
491
digitB.SetPartEventNumber(evnum);
492
digitB.SetNumberOfPes(-999);
493
digitB.SetPositionMin4(-999);
494
digitB.SetPositionMin5(-999);
495
digitB.SetPositionMin6(-999);
496
digitB.SetPositionMin7(-999);
497
digitB.SetPulse4(-999);
498
digitB.SetPulse5(-999);
499
digitB.SetPulse6(-999);
500
digitB.SetPulse7(-999);
501
digitB.SetCoincidences(-999);
502
503
digitA.SetArrivalTime0(std::accumulate(arrival_time0_each_hit_a.begin(),
504
arrival_time0_each_hit_a.end(),
505
0.0) / arrival_time0_each_hit_a.size());
506
digitA.SetArrivalTime1(std::accumulate(arrival_time1_each_hit_a.begin(),
507
arrival_time1_each_hit_a.end(),
508
0.0) / arrival_time1_each_hit_a.size());
509
digitA.SetArrivalTime2(std::accumulate(arrival_time2_each_hit_a.begin(),
510
arrival_time2_each_hit_a.end(),
511
0.0) / arrival_time2_each_hit_a.size());
512
digitA.SetArrivalTime3(std::accumulate(arrival_time3_each_hit_a.begin(),
513
arrival_time3_each_hit_a.end(),
514
0.0) / arrival_time3_each_hit_a.size());
515
double total_npe = std::accumulate(npe_each_hit_a.begin(),
516
npe_each_hit_a.end(),
517
0.0);
518
digitA.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
519
digitA.SetPartEventNumber(evnum);
520
digitA.SetNumberOfPes(total_npe);
521
digitA.SetPositionMin0(0);
522
digitA.SetPositionMin1(0);
523
digitA.SetPositionMin2(0);
524
digitA.SetPositionMin3(0);
525
digitA.SetPulse0(0);
526
digitA.SetPulse1(0);
527
digitA.SetPulse2(0);
528
digitA.SetPulse3(0);
529
digitA.SetCoincidences(0);
530
} else {
531
digitA.SetArrivalTime0(std::accumulate(arrival_time0_each_hit_a.begin(),
532
arrival_time0_each_hit_a.end(),
533
0.0) / arrival_time0_each_hit_a.size());
534
digitA.SetArrivalTime1(std::accumulate(arrival_time1_each_hit_a.begin(),
535
arrival_time1_each_hit_a.end(),
536
0.0) / arrival_time1_each_hit_a.size());
537
digitA.SetArrivalTime2(std::accumulate(arrival_time2_each_hit_a.begin(),
538
arrival_time2_each_hit_a.end(),
539
0.0) / arrival_time2_each_hit_a.size());
540
digitA.SetArrivalTime3(std::accumulate(arrival_time3_each_hit_a.begin(),
541
arrival_time3_each_hit_a.end(),
542
0.0) / arrival_time3_each_hit_a.size());
543
double total_npe = std::accumulate(npe_each_hit_a.begin(),
544
npe_each_hit_a.end(),
545
0.0);
546
digitA.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
547
digitA.SetPartEventNumber(evnum);
548
digitA.SetNumberOfPes(total_npe);
549
digitA.SetPositionMin0(0);
550
digitA.SetPositionMin1(0);
551
digitA.SetPositionMin2(0);
552
digitA.SetPositionMin3(0);
553
digitA.SetPulse0(0);
554
digitA.SetPulse1(0);
555
digitA.SetPulse2(0);
556
digitA.SetPulse3(0);
557
digitA.SetCoincidences(0);
558
digitB.SetArrivalTime4(std::accumulate(arrival_time4_each_hit_b.begin(),
559
arrival_time4_each_hit_b.end(),
560
0.0) / arrival_time4_each_hit_b.size());
561
digitB.SetArrivalTime5(std::accumulate(arrival_time5_each_hit_b.begin(),
562
arrival_time5_each_hit_b.end(),
563
0.0) / arrival_time5_each_hit_b.size());
564
digitB.SetArrivalTime6(std::accumulate(arrival_time6_each_hit_b.begin(),
565
arrival_time6_each_hit_b.end(),
566
0.0) / arrival_time6_each_hit_b.size());
567
digitB.SetArrivalTime7(std::accumulate(arrival_time7_each_hit_b.begin(),
568
arrival_time7_each_hit_b.end(),
569
0.0) / arrival_time7_each_hit_b.size());
570
total_npe = std::accumulate(npe_each_hit_b.begin(),
571
npe_each_hit_b.end(),
572
0.0);
573
digitB.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
574
digitB.SetPartEventNumber(evnum);
575
digitB.SetNumberOfPes(total_npe);
576
digitB.SetPositionMin4(0);
577
digitB.SetPositionMin5(0);
578
digitB.SetPositionMin6(0);
579
digitB.SetPositionMin7(0);
580
digitB.SetPulse4(0);
581
digitB.SetPulse5(0);
582
digitB.SetPulse6(0);
583
digitB.SetPulse7(0);
584
digitB.SetCoincidences(0);
585
}
586
587
// Now fill in the digits:
588
CkovDigit aDigit;
589
aDigit.SetCkovA(digitA);
590
aDigit.SetCkovB(digitB);
591
digit_array->push_back(aDigit);
497
592
}
498
593
//=====================================================================
499
594
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1