~ajdobbs/maus/event-viewer
» Revision
910.1.11
: src/map/MapCppCkovMCDigitizer/MapCppCkovMCDigitizer.cc
« back to all changes in this revision
Viewing changes to src/map/MapCppCkovMCDigitizer/MapCppCkovMCDigitizer.cc
- browse files at revision 910.1.11
- compare with another revision
- download diff
- download tarball
- view history from revision 910.1.11
- Committer: Durga Rajaram
- Date: 2015-09-02 11:26:55 UTC
- mfrom: (896.3.4 maus)
- mto: This revision was merged to the branch mainline in revision 921.
- Revision ID: durga@fnal.gov-20150902112655-j2pwk1wlno8afsic
merge in changes from Ao
added
removed
src/map/MapCppCkovMCDigitizer/MapCppCkovMCDigitizer.cc
17
17
18
18
#include "src/map/MapCppCkovMCDigitizer/MapCppCkovMCDigitizer.hh"
19
19
#include <string>
20
#include <numeric>
20
21
#include "Config/MiceModule.hh"
21
22
#include "Utils/Exception.hh"
22
23
#include "API/PyWrapMapBase.hh"
120
121
// loop over events
121
122
for ( unsigned int i = 0; i < mcEvts->size(); i++ ) {
122
123
// Get Ckov hits for this event
123
CkovHitArray *hits = spill->GetAnMCEvent(i).GetCkovHits();
124
CkovHitArray *hit_array = spill->GetAnMCEvent(i).GetCkovHits();
125
if (fDebug) printf("\nNumber of hits: %i\n", hit_array->size());
124
126
125
127
double gentime = spill->GetAnMCEvent(i).GetPrimary()->GetTime();
126
128
CkovEvent *evt = new CkovEvent();
127
129
(*recEvts)[i]->SetCkovEvent(evt);
128
130
129
131
// process only if there are hits
130
// if _hits.size() = 0, the make_digits and fill_ckov_evt functions will
131
// return null. The Ckov recon expects something - either null or data, can't just skip
132
if (hits) {
132
if (hit_array) {
133
133
134
134
// pick out ckov hits, digitize and store
135
CkovTmpDigits tmpAllDigits = make_ckov_digits(hits, gentime);
136
137
CkovDigitArray* CkovDigArray = evt->GetCkovDigitArrayPtr();
138
fill_ckov_evt(i, tmpAllDigits, CkovDigArray);
139
140
// DR 2015-08-31 -- the loops below are now done within fill_ckov_evt
141
// ------>NOW, go through the stations and fill Ckov events
142
// loop over ckov stations
143
/*
144
Json::Value ckov_digitAB(Json::arrayValue);
145
ckov_digitAB.clear();
146
for (int snum = 0; snum < 2; ++snum) {
147
for (unsigned int kk = 0; kk < all_ckov_digits.size(); ++kk) {
148
// check that this digit belongs to the station we are trying to fill
149
if (all_ckov_digits[kk]["station"].asInt() != snum) continue;
150
ckov_digitAB.append(fill_ckov_evt(i, snum, all_ckov_digits, kk));
151
}
152
} // end loop over stations
153
root["recon_events"][i]["ckov_event"]["ckov_digits"] = ckov_digitAB;
154
*/
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);
155
139
} // end check-if-hits
156
140
} // end loop over events
157
141
}
158
142
159
143
//////////////////////////////////////////////////////////////////////
160
CkovTmpDigits MapCppCkovMCDigitizer::make_ckov_digits(CkovHitArray* hits,
161
double gentime) const {
162
163
CkovTmpDigits tmpDigits(0);
164
if (hits->size() == 0) return tmpDigits;
165
166
for (unsigned int j = 0; j < hits->size(); ++j) { // j-th hit
167
CkovHit hit = hits->at(j);
144
multi_ckov_dig MapCppCkovMCDigitizer::make_ckov_digits(CkovHitArray* hit_array,
145
double gentime) const {
146
147
multi_ckov_dig all_ckov_dig(0);
148
if (hit_array->size() == 0) return all_ckov_dig;
149
150
for (unsigned int j = 0; j < hit_array->size(); ++j) { // j-th hit
151
CkovHit hit = hit_array->at(j);
168
152
if (fDebug) std::cout << "=================== hit# " << j << std::endl;
169
153
170
154
// make sure we can get the station number
251
235
<< " " << tdc2 << " " << tdc3 << std::endl;
252
236
}
253
237
254
aTmpCkovDigit aTmpDigit;
255
aTmpDigit.fChannelId = hit.GetChannelId();
256
aTmpDigit.fMom = hit.GetMomentum();
257
aTmpDigit.fTime = hit.GetTime();
258
aTmpDigit.fPos = hit.GetPosition();
259
aTmpDigit.fIsUsed = 0;
260
aTmpDigit.fNpe = npe;
238
one_ckov_dig a_ckov_dig;
239
a_ckov_dig.fChannelId = hit.GetChannelId();
240
a_ckov_dig.fMom = hit.GetMomentum();
241
a_ckov_dig.fTime = hit.GetTime();
242
a_ckov_dig.fPos = hit.GetPosition();
243
a_ckov_dig.fIsUsed = 0;
244
a_ckov_dig.fNpe = npe;
261
245
262
246
// set the key for this channel
263
aTmpDigit.fStation = stn;
264
aTmpDigit.fNpe0 = npe0;
265
aTmpDigit.fLeadingTime0 = tdc0;
266
aTmpDigit.fRawTime0 = time0;
267
aTmpDigit.fNpe1 = npe1;
268
aTmpDigit.fLeadingTime1 = tdc1;
269
aTmpDigit.fRawTime1 = time1;
270
aTmpDigit.fNpe2 = npe2;
271
aTmpDigit.fLeadingTime2 = tdc2;
272
aTmpDigit.fRawTime2 = time2;
273
aTmpDigit.fNpe3 = npe3;
274
aTmpDigit.fLeadingTime3 = tdc3;
275
aTmpDigit.fRawTime3 = time3;
276
tmpDigits.push_back(aTmpDigit);
277
} // for (unsigned int j = 0; j < hits.size(); ++j)
247
a_ckov_dig.fStation = stn;
248
a_ckov_dig.fNpe0 = npe0;
249
a_ckov_dig.fLeadingTime0 = tdc0;
250
a_ckov_dig.fRawTime0 = time0;
251
a_ckov_dig.fNpe1 = npe1;
252
a_ckov_dig.fLeadingTime1 = tdc1;
253
a_ckov_dig.fRawTime1 = time1;
254
a_ckov_dig.fNpe2 = npe2;
255
a_ckov_dig.fLeadingTime2 = tdc2;
256
a_ckov_dig.fRawTime2 = time2;
257
a_ckov_dig.fNpe3 = npe3;
258
a_ckov_dig.fLeadingTime3 = tdc3;
259
a_ckov_dig.fRawTime3 = time3;
260
all_ckov_dig.push_back(a_ckov_dig);
261
} // for (unsigned int j = 0; j < hit_array.size(); ++j)
278
262
279
return tmpDigits;
263
return all_ckov_dig;
280
264
}
281
265
//////////////////////////////////////////////////////////////////////
282
266
double MapCppCkovMCDigitizer::get_npe(CkovHit& hit) const {
283
267
284
268
double nphot = 0.;
285
269
double nptmp = 0.;
286
/* these are defined in the header & initialized at birth - DR 2015-08-31
287
double muon_thrhold = 0.;
288
double charge_per_pe = 0.;
289
// scaling factor to data
290
double scaling = 0.935;
291
*/
292
270
// mass of the hit particle;
293
271
double particle_mass = hit.GetMass();
294
272
// momentum of the hit particle;
299
277
300
278
int hitStation = hit.GetChannelId()->GetStation();
301
279
int hitPid = hit.GetParticleId();
302
/* the initialization of these is now done in birth so as to keep
303
* all hardcoded values in one place - DR 2015-08-31
304
if (hitStation == 0) {
305
muon_thrhold = 275.0;
306
charge_per_pe = 17.3;
307
} else if (hitStation == 1) {
308
muon_thrhold = 210.0;
309
charge_per_pe = 26;
310
}
311
*/
312
280
313
281
// momentum threshold, see Lucien's Python code.
314
282
double pp_threshold = muon_thrhold[hitStation] * particle_mass / 105.6;
407
375
408
376
//////////////////////////////////////////////////////////////////////
409
377
void MapCppCkovMCDigitizer::fill_ckov_evt(int evnum,
410
CkovTmpDigits& tmpAllDigits,
411
CkovDigitArray* CkovDigArray) const {
412
// return null if this evt had no ckov hits
413
if (tmpAllDigits.size() == 0) return;
414
415
double npe;
416
int ndigs = 0;
417
418
for (int snum = 0; snum < 2; ++snum) {
419
for (unsigned int kk = 0; kk < tmpAllDigits.size(); ++kk) {
420
// check that this hit has not already been used/filled
421
if (tmpAllDigits[kk].fStation != snum) continue;
422
CkovDigit aDigit;
423
if (tmpAllDigits[kk].fIsUsed == 0) {
424
ndigs++;
425
npe = tmpAllDigits[kk].fNpe;
426
// 23 is the ADC factor.
427
// the factor - ckovNpeToAdc is set at birth, defined in header
428
// keeping for-now-hardcoded definitions in one place
429
// -- DR 2015-08-31
430
int adc = static_cast<int>(npe * ckovNpeToAdc);
431
432
CkovA digitA;
433
CkovB digitB;
434
if (snum == 0) {
435
digitA.SetArrivalTime0(tmpAllDigits[kk].fLeadingTime0);
436
digitA.SetArrivalTime1(tmpAllDigits[kk].fLeadingTime1);
437
digitA.SetArrivalTime2(tmpAllDigits[kk].fLeadingTime2);
438
digitA.SetArrivalTime3(tmpAllDigits[kk].fLeadingTime3);
439
digitA.SetTotalCharge(adc);
440
digitA.SetPartEventNumber(evnum);
441
digitA.SetNumberOfPes(npe);
442
443
// Does not matter too much but..
444
// I think the pulse0/1/2/3 should be adc/4 -- DR 2015-08-31
445
digitA.SetPositionMin0(0);
446
digitA.SetPositionMin1(0);
447
digitA.SetPositionMin2(0);
448
digitA.SetPositionMin3(0);
449
digitA.SetPulse0(0);
450
digitA.SetPulse1(0);
451
digitA.SetPulse2(0);
452
digitA.SetPulse3(0);
453
digitA.SetCoincidences(0);
454
aDigit.SetCkovA(digitA);
455
} else {
456
digitB.SetArrivalTime4(tmpAllDigits[kk].fLeadingTime0);
457
digitB.SetArrivalTime5(tmpAllDigits[kk].fLeadingTime1);
458
digitB.SetArrivalTime6(tmpAllDigits[kk].fLeadingTime2);
459
digitB.SetArrivalTime7(tmpAllDigits[kk].fLeadingTime3);
460
digitB.SetTotalCharge(adc);
461
digitB.SetPartEventNumber(evnum);
462
digitB.SetNumberOfPes(npe);
463
464
// Does not matter too much but..
465
// I think the pulse4/5/6/7 should be adc/4 -- DR 2015-08-31
466
digitB.SetPositionMin4(0);
467
digitB.SetPositionMin5(0);
468
digitB.SetPositionMin6(0);
469
digitB.SetPositionMin7(0);
470
digitB.SetPulse4(0);
471
digitB.SetPulse5(0);
472
digitB.SetPulse6(0);
473
digitB.SetPulse7(0);
474
digitB.SetCoincidences(0);
475
aDigit.SetCkovB(digitB);
476
}
477
tmpAllDigits[kk].fIsUsed = true;
478
} // end check if-digit-used
479
CkovDigArray->push_back(aDigit);
480
} // end loop over tmp digits
481
} // end loop over stations
378
multi_ckov_dig& all_ckov_dig,
379
CkovDigitArray* digit_array) const {
380
// return null if this evt had no ckov hits, no need to go further;
381
if (all_ckov_dig.size() == 0) return;
382
383
/* If there are hits:
384
* For this given event number, there might be multiple hits;
385
* If the hits belong to 2 stations, they should be filled simultaneously to
386
* the digits, but in their respective structure;
387
* Otherwise, we need to integrate the digits together and get an average;
388
*/
389
std::vector<double> npe_each_hit_a;
390
std::vector<double> npe_each_hit_b;
391
std::vector<double> arrival_time0_each_hit_a;
392
std::vector<double> arrival_time1_each_hit_a;
393
std::vector<double> arrival_time2_each_hit_a;
394
std::vector<double> arrival_time3_each_hit_a;
395
std::vector<double> arrival_time4_each_hit_b;
396
std::vector<double> arrival_time5_each_hit_b;
397
std::vector<double> arrival_time6_each_hit_b;
398
std::vector<double> arrival_time7_each_hit_b;
399
400
// Count how many hits each station has
401
int num_of_hits_a = 0;
402
int num_of_hits_b = 0;
403
404
// Station number
405
int 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
snum = all_ckov_dig[kk].fStation;
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
434
// Now, if either num_of_hits_a(b) == 0, write the default values to the station:
435
// Can not be 0 at the same time. If there are no hits, it won't reach this far.
436
CkovA digitA;
437
CkovB digitB;
438
if (num_of_hits_a == 0) {
439
440
digitA.SetArrivalTime0(-999);
441
digitA.SetArrivalTime1(-999);
442
digitA.SetArrivalTime2(-999);
443
digitA.SetArrivalTime3(-999);
444
digitA.SetTotalCharge(-999);
445
digitA.SetPartEventNumber(evnum);
446
digitA.SetNumberOfPes(-999);
447
digitA.SetPositionMin0(-999);
448
digitA.SetPositionMin1(-999);
449
digitA.SetPositionMin2(-999);
450
digitA.SetPositionMin3(-999);
451
digitA.SetPulse0(-999);
452
digitA.SetPulse1(-999);
453
digitA.SetPulse2(-999);
454
digitA.SetPulse3(-999);
455
digitA.SetCoincidences(-999);
456
457
digitB.SetArrivalTime4(std::accumulate(arrival_time4_each_hit_b.begin(),
458
arrival_time4_each_hit_b.end(),
459
0.0) / arrival_time4_each_hit_b.size());
460
digitB.SetArrivalTime5(std::accumulate(arrival_time5_each_hit_b.begin(),
461
arrival_time5_each_hit_b.end(),
462
0.0) / arrival_time5_each_hit_b.size());
463
digitB.SetArrivalTime6(std::accumulate(arrival_time6_each_hit_b.begin(),
464
arrival_time6_each_hit_b.end(),
465
0.0) / arrival_time6_each_hit_b.size());
466
digitB.SetArrivalTime7(std::accumulate(arrival_time7_each_hit_b.begin(),
467
arrival_time7_each_hit_b.end(),
468
0.0) / arrival_time7_each_hit_b.size());
469
double total_npe = std::accumulate(npe_each_hit_b.begin(),
470
npe_each_hit_b.end(),
471
0.0);
472
digitB.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
473
digitB.SetPartEventNumber(evnum);
474
digitB.SetNumberOfPes(total_npe);
475
digitB.SetPositionMin4(0);
476
digitB.SetPositionMin5(0);
477
digitB.SetPositionMin6(0);
478
digitB.SetPositionMin7(0);
479
digitB.SetPulse4(0);
480
digitB.SetPulse5(0);
481
digitB.SetPulse6(0);
482
digitB.SetPulse7(0);
483
digitB.SetCoincidences(0);
484
} else if (num_of_hits_b == 0) {
485
digitB.SetArrivalTime4(-999);
486
digitB.SetArrivalTime5(-999);
487
digitB.SetArrivalTime6(-999);
488
digitB.SetArrivalTime7(-999);
489
digitB.SetTotalCharge(-999);
490
digitB.SetPartEventNumber(evnum);
491
digitB.SetNumberOfPes(-999);
492
digitB.SetPositionMin4(-999);
493
digitB.SetPositionMin5(-999);
494
digitB.SetPositionMin6(-999);
495
digitB.SetPositionMin7(-999);
496
digitB.SetPulse4(-999);
497
digitB.SetPulse5(-999);
498
digitB.SetPulse6(-999);
499
digitB.SetPulse7(-999);
500
digitB.SetCoincidences(-999);
501
digitA.SetArrivalTime0(std::accumulate(arrival_time0_each_hit_a.begin(),
502
arrival_time0_each_hit_a.end(),
503
0.0) / arrival_time0_each_hit_a.size());
504
digitA.SetArrivalTime1(std::accumulate(arrival_time1_each_hit_a.begin(),
505
arrival_time1_each_hit_a.end(),
506
0.0) / arrival_time1_each_hit_a.size());
507
digitA.SetArrivalTime2(std::accumulate(arrival_time2_each_hit_a.begin(),
508
arrival_time2_each_hit_a.end(),
509
0.0) / arrival_time2_each_hit_a.size());
510
digitA.SetArrivalTime3(std::accumulate(arrival_time3_each_hit_a.begin(),
511
arrival_time3_each_hit_a.end(),
512
0.0) / arrival_time3_each_hit_a.size());
513
double total_npe = std::accumulate(npe_each_hit_a.begin(),
514
npe_each_hit_a.end(),
515
0.0);
516
digitA.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
517
digitA.SetPartEventNumber(evnum);
518
digitA.SetNumberOfPes(total_npe);
519
digitA.SetPositionMin0(0);
520
digitA.SetPositionMin1(0);
521
digitA.SetPositionMin2(0);
522
digitA.SetPositionMin3(0);
523
digitA.SetPulse0(0);
524
digitA.SetPulse1(0);
525
digitA.SetPulse2(0);
526
digitA.SetPulse3(0);
527
digitA.SetCoincidences(0);
528
} else {
529
digitA.SetArrivalTime0(std::accumulate(arrival_time0_each_hit_a.begin(),
530
arrival_time0_each_hit_a.end(),
531
0.0) / arrival_time0_each_hit_a.size());
532
digitA.SetArrivalTime1(std::accumulate(arrival_time1_each_hit_a.begin(),
533
arrival_time1_each_hit_a.end(),
534
0.0) / arrival_time1_each_hit_a.size());
535
digitA.SetArrivalTime2(std::accumulate(arrival_time2_each_hit_a.begin(),
536
arrival_time2_each_hit_a.end(),
537
0.0) / arrival_time2_each_hit_a.size());
538
digitA.SetArrivalTime3(std::accumulate(arrival_time3_each_hit_a.begin(),
539
arrival_time3_each_hit_a.end(),
540
0.0) / arrival_time3_each_hit_a.size());
541
double total_npe = std::accumulate(npe_each_hit_a.begin(),
542
npe_each_hit_a.end(),
543
0.0);
544
digitA.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
545
digitA.SetPartEventNumber(evnum);
546
digitA.SetNumberOfPes(total_npe);
547
digitA.SetPositionMin0(0);
548
digitA.SetPositionMin1(0);
549
digitA.SetPositionMin2(0);
550
digitA.SetPositionMin3(0);
551
digitA.SetPulse0(0);
552
digitA.SetPulse1(0);
553
digitA.SetPulse2(0);
554
digitA.SetPulse3(0);
555
digitA.SetCoincidences(0);
556
digitB.SetArrivalTime4(std::accumulate(arrival_time4_each_hit_b.begin(),
557
arrival_time4_each_hit_b.end(),
558
0.0) / arrival_time4_each_hit_b.size());
559
digitB.SetArrivalTime5(std::accumulate(arrival_time5_each_hit_b.begin(),
560
arrival_time5_each_hit_b.end(),
561
0.0) / arrival_time5_each_hit_b.size());
562
digitB.SetArrivalTime6(std::accumulate(arrival_time6_each_hit_b.begin(),
563
arrival_time6_each_hit_b.end(),
564
0.0) / arrival_time6_each_hit_b.size());
565
digitB.SetArrivalTime7(std::accumulate(arrival_time7_each_hit_b.begin(),
566
arrival_time7_each_hit_b.end(),
567
0.0) / arrival_time7_each_hit_b.size());
568
total_npe = std::accumulate(npe_each_hit_b.begin(),
569
npe_each_hit_b.end(),
570
0.0);
571
digitB.SetTotalCharge(static_cast<int> (total_npe * ckovNpeToAdc));
572
digitB.SetPartEventNumber(evnum);
573
digitB.SetNumberOfPes(total_npe);
574
digitB.SetPositionMin4(0);
575
digitB.SetPositionMin5(0);
576
digitB.SetPositionMin6(0);
577
digitB.SetPositionMin7(0);
578
digitB.SetPulse4(0);
579
digitB.SetPulse5(0);
580
digitB.SetPulse6(0);
581
digitB.SetPulse7(0);
582
digitB.SetCoincidences(0);
583
}
584
585
// Now fill in the digits:
586
CkovDigit aDigit;
587
aDigit.SetCkovA(digitA);
588
aDigit.SetCkovB(digitB);
589
digit_array->push_back(aDigit);
482
590
}
483
591
//=====================================================================
484
592
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1