~durga/maus/rel709

« back to all changes in this revision

Viewing changes to src/map/MapCppEMRMCDigitization/MapCppEMRMCDigitization.cc

  • Committer: Adam Dobbs
  • Date: 2014-12-11 13:26:05 UTC
  • mfrom: (659.1.96 release-candidate)
  • Revision ID: phuccj@gmail.com-20141211132605-6g6j2927elggi2q6
Tags: MAUS-v0.9.2
MAUS-v0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
 
2
 *
 
3
 * MAUS is free software: you can redistribute it and/or modify
 
4
 * it under the terms of the GNU General Public License as published by
 
5
 * the Free Software Foundation, either version 3 of the License, or
 
6
 * (at your option) any later version.
 
7
 *
 
8
 * MAUS is distributed in the hope that it will be useful,
 
9
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
 * GNU General Public License for more details.
 
12
 *
 
13
 * You should have received a copy of the GNU General Public License
 
14
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
 
15
 *
 
16
 */
 
17
 
 
18
#include "Utils/CppErrorHandler.hh"
 
19
#include "Utils/JsonWrapper.hh"
 
20
#include "Interface/dataCards.hh"
 
21
#include "DataStructure/MCEvent.hh"
 
22
#include "DataStructure/Spill.hh"
 
23
#include "API/PyWrapMapBase.hh"
 
24
#include "Converter/DataConverters/CppJsonSpillConverter.hh"
 
25
#include "Converter/DataConverters/JsonCppSpillConverter.hh"
 
26
 
 
27
#include "src/map/MapCppEMRMCDigitization/MapCppEMRMCDigitization.hh"
 
28
 
 
29
namespace MAUS {
 
30
 
 
31
PyMODINIT_FUNC init_MapCppEMRMCDigitization(void) {
 
32
  PyWrapMapBase<MAUS::MapCppEMRMCDigitization>::PyWrapMapBaseModInit
 
33
                                                ("MapCppEMRMCDigitization", "", "", "", "");
 
34
}
 
35
 
 
36
MapCppEMRMCDigitization::MapCppEMRMCDigitization()
 
37
    : MapBase<MAUS::Data>("MapCppEMRMCDigitization") {
 
38
}
 
39
 
 
40
void MapCppEMRMCDigitization::_birth(const std::string& argJsonConfigDocument) {
 
41
  _classname = "MapCppEMRMCDigitization";
 
42
  char* pMAUS_ROOT_DIR = getenv("MAUS_ROOT_DIR");
 
43
 
 
44
  if (!pMAUS_ROOT_DIR) {
 
45
    throw MAUS::Exception(Exception::recoverable,
 
46
                      "Could not resolve ${MAUS_ROOT_DIR} environment variable",
 
47
                      "MapCppEMRMCDigitization::birth");
 
48
  }
 
49
 
 
50
  //  JsonCpp setup
 
51
  Json::Value configJSON;
 
52
  configJSON = JsonWrapper::StringToJson(argJsonConfigDocument);
 
53
 
 
54
  // Fetch variables
 
55
  _number_of_planes = configJSON["EMRnumberOfPlanes"].asInt();
 
56
  _number_of_bars = configJSON["EMRnumberOfBars"].asInt();
 
57
 
 
58
  _tot_noise_up = configJSON["EMRtotNoiseUp"].asInt();
 
59
  _tot_noise_low = configJSON["EMRtotNoiseLow"].asInt();
 
60
 
 
61
  _do_sampling = configJSON["EMRdoSampling"].asInt();
 
62
  _nph_per_MeV = configJSON["EMRnphPerMeV"].asInt();
 
63
  _seed = configJSON["EMRseed"].asInt();
 
64
  _trap_eff = configJSON["EMRtrapEff"].asDouble();
 
65
  _QE_MAPMT = configJSON["EMRqeMAPMT"].asDouble();
 
66
  _QE_SAPMT = configJSON["EMRqeSAPMT"].asDouble();
 
67
  _nADC_per_pe_MAPMT = configJSON["EMRnadcPerPeMAPMT"].asDouble();
 
68
  _nADC_per_pe_SAPMT = configJSON["EMRnadcPerPeSAPMT"].asDouble();
 
69
  _electronics_response_spread_MAPMT = configJSON["EMRelectronicsResponseSpreadMAPMT"].asDouble();
 
70
  _electronics_response_spread_SAPMT = configJSON["EMRelectronicsResponseSpreadSAPMT"].asDouble();
 
71
  _dbb_count = configJSON["EMRdbbCount"].asDouble();
 
72
  _fadc_count = configJSON["EMRfadcCount"].asDouble();
 
73
  _time_response_spread = configJSON["EMRtimeResponseSpread"].asInt();
 
74
  _tot_func_p1 = configJSON["EMRtotFuncP1"].asDouble();
 
75
  _tot_func_p2 = configJSON["EMRtotFuncP2"].asDouble();
 
76
  _tot_func_p3 = configJSON["EMRtotFuncP3"].asDouble();
 
77
  _tot_func_p4 = configJSON["EMRtotFuncP4"].asDouble();
 
78
  _acquisition_window = configJSON["EMRacquisitionWindow"].asInt();
 
79
  _signal_integration_window = configJSON["EMRsignalIntegrationWindow"].asInt();
 
80
  _arrival_time_shift = configJSON["EMRarrivalTimeShift"].asInt();
 
81
  _arrival_time_gaus_width = configJSON["EMRarrivalTimeGausWidth"].asDouble();
 
82
  _arrival_time_uniform_width = configJSON["EMRarrivalTimeUniformWidth"].asDouble();
 
83
  _pulse_shape_landau_width = configJSON["EMRpulseShapeLandauWidth"].asDouble();
 
84
  _fom = configJSON["EMRfom"].asString();
 
85
  _birks_constant = configJSON["EMRbirksConstant"].asDouble();
 
86
  _average_path_length = configJSON["EMRaveragePathLength"].asDouble();
 
87
  _signal_energy_threshold = configJSON["EMRsignalEnergyThreshold"].asDouble();
 
88
  _baseline_spread = configJSON["EMRbaselineSpread"].asInt();
 
89
  _maximum_noise_level = configJSON["EMRmaximumNoiseLevel"].asInt();
 
90
  _deltat_shift = configJSON["EMRdeltatShift"].asDouble();
 
91
  _baseline_position = configJSON["EMRbaselinePosition"].asDouble();
 
92
  _arrival_time_spread = configJSON["EMRarrivalTimeSpread"].asDouble();
 
93
 
 
94
  // Generate random noise in each SAPMT
 
95
  _rand = new TRandom3(_seed);
 
96
 
 
97
  for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
98
    _baseline[iPlane] = _baseline_position
 
99
                 + static_cast<int>(_rand->Uniform(-_baseline_spread, _baseline_spread));
 
100
    _noise_level[iPlane] = static_cast<int>(_rand->Uniform(0, _maximum_noise_level));
 
101
    _noise_position[iPlane] = static_cast<int>(_rand->Uniform(0, 1.9999999));
 
102
  }
 
103
 
 
104
  // Load the EMR calibration map
 
105
  bool loaded = _calibMap.InitializeFromCards(configJSON);
 
106
  if (!loaded)
 
107
    throw(Exception(Exception::recoverable,
 
108
          "Could not find EMR calibration map",
 
109
          "MapCppEMRMCDigitizer::birth"));
 
110
 
 
111
  // Load the EMR attenuation map
 
112
  loaded = _attenMap.InitializeFromCards(configJSON);
 
113
  if (!loaded)
 
114
    throw(Exception(Exception::recoverable,
 
115
          "Could not find EMR attenuation map",
 
116
          "MapCppEMRMCDigitizer::birth"));
 
117
}
 
118
 
 
119
void MapCppEMRMCDigitization::_death() {
 
120
}
 
121
 
 
122
void MapCppEMRMCDigitization::_process(Data *data) const {
 
123
 
 
124
  // Get spill, break if there's no DAQ data
 
125
  Spill *spill = data->GetSpill();
 
126
 
 
127
  if (spill->GetMCEvents() == NULL)
 
128
      return;
 
129
 
 
130
  int nPartEvents = spill->GetMCEventSize();
 
131
 
 
132
  // Check the Recon event array is initialised, and if not make it so
 
133
  if (!spill->GetReconEvents()) {
 
134
    ReconEventPArray* revts = new ReconEventPArray();
 
135
    spill->SetReconEvents(revts);
 
136
  }
 
137
 
 
138
  // Set primary hits/noise time windows
 
139
  std::vector<double> delta_t_array;
 
140
  delta_t_array.resize(0);
 
141
 
 
142
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
143
    EMRHitArray *hits = spill->GetAnMCEvent(iPe).GetEMRHits();
 
144
 
 
145
    if (hits) {
 
146
      int nHits = hits->size();
 
147
 
 
148
      Primary *primary = spill->GetAnMCEvent(iPe).GetPrimary();
 
149
      double pTime = primary->GetTime(); // ns
 
150
 
 
151
      for (int iHit = 0; iHit < nHits; iHit++) {
 
152
        EMRHit hit = hits->at(iHit);
 
153
 
 
154
        int time = static_cast<int>(hit.GetTime()); // ns
 
155
        double delta_t = time - pTime;
 
156
 
 
157
        delta_t_array.push_back(delta_t);
 
158
      }
 
159
    }
 
160
  }
 
161
 
 
162
  if (!delta_t_array.size()) return;
 
163
 
 
164
  int lt = static_cast<int>(*std::min_element(delta_t_array.begin(), delta_t_array.end())
 
165
                            /_dbb_count);
 
166
  int deltat_limits[4] = {lt, lt+20, lt+20, lt+40};
 
167
 
 
168
  // Create the DBB and fADC arrays with n events (1 per trigger)
 
169
  EMRDBBEventVector emr_dbb_events_tmp = get_dbb_data_tmp(nPartEvents);
 
170
  EMRfADCEventVector emr_fadc_events_tmp = get_fadc_data_tmp(nPartEvents);
 
171
 
 
172
  // Fill the Recon event array with G4 Spill information
 
173
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
174
    EMRHitArray *EMRhits = spill->GetAnMCEvent(iPe).GetEMRHits();
 
175
 
 
176
    if (EMRhits) {
 
177
      Primary *primary = spill->GetAnMCEvent(iPe).GetPrimary();
 
178
      double pTime = primary->GetTime(); // ns
 
179
 
 
180
      processDBB(EMRhits, iPe, pTime, emr_dbb_events_tmp, emr_fadc_events_tmp);
 
181
      processFADC(EMRhits, iPe, emr_fadc_events_tmp);
 
182
    }
 
183
  }
 
184
 
 
185
  fill(spill, nPartEvents, emr_dbb_events_tmp, emr_fadc_events_tmp);
 
186
 
 
187
  // Reset/Resize DBB and fADC arrays with n+2 events (1 per trigger + noise + decays)
 
188
  emr_dbb_events_tmp = get_dbb_data_tmp(nPartEvents+2);
 
189
  emr_fadc_events_tmp = get_fadc_data_tmp(nPartEvents+2);
 
190
 
 
191
  // Digitize the Recon event array
 
192
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
193
    EMREvent *EMRevt = spill->GetReconEvents()->at(iPe)->GetEMREvent();
 
194
 
 
195
    if (EMRevt) {
 
196
      digitize(EMRevt, iPe, nPartEvents, deltat_limits, emr_dbb_events_tmp, emr_fadc_events_tmp);
 
197
    }
 
198
  }
 
199
 
 
200
  fill(spill, nPartEvents+2, emr_dbb_events_tmp, emr_fadc_events_tmp);
 
201
}
 
202
 
 
203
void MapCppEMRMCDigitization::processDBB(MAUS::EMRHitArray *EMRhits,
 
204
                                         int xPe,
 
205
                                         double pTime,
 
206
                                         EMRDBBEventVector& emr_dbb_events_tmp,
 
207
                                         EMRfADCEventVector& emr_fadc_events_tmp) const {
 
208
 
 
209
  int nHits = EMRhits->size();
 
210
 
 
211
  for (int iHit = 0; iHit < nHits; iHit++) {
 
212
 
 
213
    EMRHit hit = EMRhits->at(iHit);
 
214
 
 
215
    // Channel ID
 
216
    EMRChannelId *id = hit.GetChannelId();
 
217
 
 
218
    int g4barid = id->GetBar();
 
219
    int planeid = g4barid / 59;
 
220
    int barid = g4barid % 59 + 1;
 
221
 
 
222
    // Energy Deposition
 
223
    double edep = hit.GetEnergyDeposited(); /*MeV*/
 
224
 
 
225
    // Position
 
226
    ThreeVector pos = hit.GetPosition();
 
227
    double pozx = pos.x(); /*mm*/
 
228
    double pozy = pos.y(); /*mm*/
 
229
    double pozz = pos.z(); /*mm*/
 
230
 
 
231
    // Time
 
232
    double time = hit.GetTime(); /*ns*/
 
233
 
 
234
    int lt  = static_cast<int>(time);/*ns*/
 
235
    int tot = static_cast<int>(edep*1000000);/*eV*/
 
236
 
 
237
    int xPlane = planeid;
 
238
    int xBar = barid;
 
239
    int delta_t = lt-pTime;
 
240
 
 
241
    // Fill EMRDBBEvent
 
242
    EMRBarHit bHit;
 
243
    bHit.SetX(pozx);
 
244
    bHit.SetY(pozy);
 
245
    bHit.SetZ(pozz);
 
246
    bHit.SetTot(tot);
 
247
    bHit.SetHitTime(lt);
 
248
    bHit.SetDeltaT(delta_t);
 
249
    emr_dbb_events_tmp[xPe][xPlane][xBar].push_back(bHit);
 
250
  }
 
251
}
 
252
 
 
253
void MapCppEMRMCDigitization::processFADC(MAUS::EMRHitArray *EMRhits,
 
254
                                          int xPe,
 
255
                                          EMRfADCEventVector& emr_fadc_events_tmp) const {
 
256
 
 
257
  int nHits = EMRhits->size();
 
258
 
 
259
  for (int iHit = 0; iHit < nHits; iHit++) {
 
260
 
 
261
    EMRHit hit = EMRhits->at(iHit);
 
262
 
 
263
    // Plane ID
 
264
    EMRChannelId *id = hit.GetChannelId();
 
265
 
 
266
    int g4barid = id->GetBar();
 
267
    int planeid = g4barid / 59;
 
268
 
 
269
    // Energy Deposition
 
270
    double edep = hit.GetEnergyDeposited(); /*MeV*/
 
271
 
 
272
    double xCharge = edep*1000000;/*eV*/
 
273
    int xPlane = planeid;
 
274
 
 
275
    xCharge += emr_fadc_events_tmp[xPe][xPlane]._charge;
 
276
    fADCdata data;
 
277
    data._orientation = planeid % 2;
 
278
    data._charge = xCharge;
 
279
    data._time = 0;
 
280
    emr_fadc_events_tmp[xPe][xPlane] = data;
 
281
  }
 
282
}
 
283
 
 
284
void MapCppEMRMCDigitization::digitize(MAUS::EMREvent *EMRevt,
 
285
                                       int xPe,
 
286
                                       int nPartEvents,
 
287
                                       int *deltat_limits,
 
288
                                       EMRDBBEventVector& emr_dbb_events_tmp,
 
289
                                       EMRfADCEventVector& emr_fadc_events_tmp) const {
 
290
 
 
291
  int nPlHits = EMRevt->GetEMRPlaneHitArray().size();
 
292
 
 
293
  for (int iPlane = 0; iPlane < nPlHits; iPlane++) {
 
294
    EMRPlaneHit *plHit = EMRevt->GetEMRPlaneHitArray().at(iPlane);
 
295
    int xPlane = plHit->GetPlane();
 
296
 
 
297
    // Find PMT gain correction coefficient
 
298
    EMRChannelKey xAverageKey(xPlane, xPlane%2, -1, "emr");
 
299
    double epsilon_MA = _calibMap.Eps(xAverageKey, "MA");
 
300
    double epsilon_SA = _calibMap.Eps(xAverageKey, "SA");
 
301
 
 
302
    int nBars = plHit->GetEMRBarArray().size();
 
303
    int naph_SAPMT = 0;
 
304
 
 
305
    for (int iBar = 0; iBar < nBars; iBar++) {
 
306
      EMRBar *bar = plHit->GetEMRBarArray().at(iBar);
 
307
      int xBar = bar->GetBar();
 
308
      int nBarHits = bar->GetEMRBarHitArray().size();
 
309
      bool matched = false;
 
310
 
 
311
      for (int iBarHit = 0; iBarHit < nBarHits; iBarHit++) {
 
312
        EMRBarHit barHit = bar->GetEMRBarHitArray().at(iBarHit);
 
313
        EMRChannelKey xKey(xPlane, xPlane%2, xBar, "emr");
 
314
        double x  = barHit.GetX();/*mm*/
 
315
        double y  = barHit.GetY();/*mm*/
 
316
        int xTot  = barHit.GetTot();/*eV*/
 
317
        int xDeltaT = barHit.GetDeltaT();/*ns*/
 
318
        int xHitTime = barHit.GetHitTime();/*ns*/
 
319
 
 
320
        /*--------- MAPMT signal simulation -------*/
 
321
 
 
322
        // simulate electronics noise and PMT dark current
 
323
        // !!! TO DO !!!
 
324
 
 
325
        // set g4 information
 
326
        double energy = xTot/1000000.0;/*MeV*/
 
327
        int time = xDeltaT;/*ns*/
 
328
 
 
329
        if (energy < _signal_energy_threshold) continue;
 
330
 
 
331
        // convert energy to the number of scintillation photons (nsph) according to Birks's law
 
332
        // !!! TO DO !!! calculate path length
 
333
        int nsph = static_cast<int>(static_cast<double>(_nph_per_MeV*energy)  // MeV
 
334
                 /(1.0+_birks_constant  // mm/MeV
 
335
                 /_average_path_length*energy));  // MeV
 
336
 
 
337
        // sample nsph with Poisson
 
338
        if (_do_sampling) nsph = _rand->Poisson(nsph);
 
339
 
 
340
        // covert nsph to the number of trapped photons (ntph)
 
341
        int ntph = static_cast<int>(static_cast<double>(nsph*_trap_eff));
 
342
 
 
343
        // sample ntph with Poisson
 
344
        if (_do_sampling) ntph = _rand->Poisson(ntph);
 
345
 
 
346
        // reduce ntph according to fibre length (naph)
 
347
        // for MAPMT
 
348
        int naph_MAPMT = static_cast<int>(static_cast<double>(ntph)
 
349
                       /2.0*_attenMap.fibreAtten(xKey, x, y, "MA"));
 
350
 
 
351
        // for SAPMT
 
352
        int naph_SAPMT_hit = static_cast<int>(static_cast<double>(ntph)
 
353
                           /2.0*_attenMap.fibreAtten(xKey, x, y, "SA"));
 
354
 
 
355
        // apply channel attenuation map
 
356
        // for MAPMT
 
357
        naph_MAPMT = static_cast<int>(static_cast<double>(naph_MAPMT)
 
358
                   *_attenMap.connectorAtten(xKey, "MA"));
 
359
 
 
360
        // for SAPMT
 
361
        naph_SAPMT_hit = static_cast<int>(static_cast<double>(naph_SAPMT_hit)
 
362
                       *_attenMap.connectorAtten(xKey, "SA"));
 
363
 
 
364
        if (time > deltat_limits[0] && time < deltat_limits[1]) naph_SAPMT += naph_SAPMT_hit;
 
365
 
 
366
        // sample naph with Poisson
 
367
        if (_do_sampling) naph_MAPMT = _rand->Poisson(naph_MAPMT);
 
368
 
 
369
        // simulate cross-talk and misalignment
 
370
        // !!! TO DO !!!
 
371
 
 
372
        // convert naph to the number of photoelectrons (npe)
 
373
        int npe = static_cast<int>(static_cast<double>(naph_MAPMT)*_QE_MAPMT);
 
374
 
 
375
        // sample npe with Poisson
 
376
        if (_do_sampling) npe = _rand->Poisson(npe);
 
377
 
 
378
        // correct npe for the photocathode non-uniformity
 
379
        // !!! TO DO !!!
 
380
 
 
381
        // correct npe for pmt gain difference
 
382
        npe = static_cast<int>(static_cast<double>(npe)*epsilon_MA);
 
383
 
 
384
        if (npe == 0) continue;
 
385
 
 
386
        // convert npe to the number of ADC counts
 
387
        int nADC = static_cast<int>(static_cast<double>(npe)*_nADC_per_pe_MAPMT);
 
388
 
 
389
        if (nADC == 0) continue;
 
390
 
 
391
        // simulate electronics response
 
392
        if (_do_sampling) nADC = static_cast<int>(_rand->Gaus(nADC, _electronics_response_spread_MAPMT));
 
393
 
 
394
        // convert nADC to Tot
 
395
        int xTotDigi = static_cast<int>(_tot_func_p1+_tot_func_p2
 
396
                     *log(static_cast<double>(nADC)
 
397
                     /_tot_func_p4+_tot_func_p3));
 
398
 
 
399
        if (xTotDigi <= 0) continue;
 
400
 
 
401
        // convert g4 time to deltaT time in ADC counts
 
402
        int xDeltaTDigi = static_cast<int>(static_cast<double>(time)/_dbb_count);
 
403
 
 
404
        // sample deltaT
 
405
        if (_do_sampling) xDeltaTDigi = static_cast<int>(_rand->Gaus(xDeltaTDigi,
 
406
                                                                     _time_response_spread));
 
407
        if (xDeltaTDigi < 0) xDeltaTDigi = 0;
 
408
 
 
409
        // correct deltaT with fibre length
 
410
        xDeltaTDigi += static_cast<int>(static_cast<double>(_attenMap.fibreDelay(xKey, x, y, "MA"))
 
411
                                        /_dbb_count);
 
412
 
 
413
        // set hit time
 
414
        int xHitTimeDigi = static_cast<int>(static_cast<double>(xHitTime)/_dbb_count);
 
415
 
 
416
        // set bar hit
 
417
        EMRBarHit bHit;
 
418
        bHit.SetTot(xTotDigi);
 
419
        bHit.SetDeltaT(xDeltaTDigi);
 
420
        bHit.SetHitTime(xHitTimeDigi);
 
421
 
 
422
        // discriminate noise and decays from signal events
 
423
        if (xDeltaTDigi >= deltat_limits[0] && xDeltaTDigi < deltat_limits[1]) {
 
424
          emr_dbb_events_tmp[xPe][xPlane][xBar].push_back(bHit);
 
425
          matched = true;
 
426
        } else if (xDeltaTDigi >= deltat_limits[2] && xDeltaTDigi < deltat_limits[3] &&
 
427
                   xTotDigi > _tot_noise_low && xTotDigi < _tot_noise_up ) {
 
428
          emr_dbb_events_tmp[nPartEvents][xPlane][xBar].push_back(bHit);
 
429
          matched = true;
 
430
        } else if (iBarHit == nBarHits-1 && !matched) {
 
431
          bHit.SetDeltaT(0);
 
432
          emr_dbb_events_tmp[nPartEvents+1][xPlane][xBar].push_back(bHit);
 
433
        }
 
434
      }
 
435
    }
 
436
 
 
437
    int xOri = plHit->GetOrientation();
 
438
 
 
439
    /*------------ SAPMT signal simulation -------*/
 
440
 
 
441
    // simulate PMT dark current
 
442
    // !!! TO DO !!!
 
443
 
 
444
    // sample naph with Poisson
 
445
    if (_do_sampling) naph_SAPMT = _rand->Poisson(naph_SAPMT);
 
446
 
 
447
    // convert naph to the number of photoelectrons (npe)
 
448
    int npe = static_cast<int>(static_cast<double>(naph_SAPMT)*_QE_SAPMT);
 
449
 
 
450
    // sample npe with Poisson
 
451
    if (_do_sampling) npe = _rand->Poisson(npe);
 
452
 
 
453
    // correct npe for pmt gain difference
 
454
    npe = static_cast<int>(static_cast<double>(npe)*epsilon_SA);
 
455
 
 
456
    // convert npe to the number of ADC counts
 
457
    int nADC = static_cast<int>(static_cast<double>(npe)*_nADC_per_pe_SAPMT);
 
458
 
 
459
    // simulate electronics response
 
460
    if (_do_sampling) nADC = static_cast<int>(_rand->Gaus(nADC,
 
461
                                                          _electronics_response_spread_SAPMT));
 
462
 
 
463
    // correct non physical values
 
464
    if (nADC < 0) nADC = 0;
 
465
 
 
466
    // set signal baseline (8bit ADC)
 
467
    int baseline = _baseline[xPlane];
 
468
 
 
469
    // set electronics noise level - number of fluctuations within acquisition window
 
470
    int noise_level = _noise_level[xPlane];
 
471
 
 
472
    // set electronics noise position (upwards/downwards fluctuations)
 
473
    int noise_position = _noise_position[xPlane];
 
474
 
 
475
    // simulate negative voltage pulse with random electronics noise
 
476
    TH1F *pulse_shape = new TH1F("pulse_shape", "pulse_shape",
 
477
                                 _acquisition_window, 0, _acquisition_window);
 
478
 
 
479
    double arrival_time = _signal_integration_window*2 + _arrival_time_shift;
 
480
 
 
481
    if (static_cast<int>(_rand->Uniform(0, 1.9999999))) arrival_time += _arrival_time_spread;
 
482
 
 
483
    arrival_time += _rand->Gaus(0, _arrival_time_gaus_width)
 
484
                 +_rand->Uniform(-_arrival_time_uniform_width, _arrival_time_uniform_width);
 
485
 
 
486
    for (int i = 0; i < noise_level; i++) {
 
487
      double noise = _rand->Uniform(0, _acquisition_window);
 
488
      int bin = static_cast<Int_t>(noise+1);
 
489
      if (pulse_shape->GetBinContent(bin) != 1) pulse_shape->Fill(noise);
 
490
    }
 
491
 
 
492
    if (noise_position) for (Int_t i = 1; i <= _acquisition_window; i++) {
 
493
      pulse_shape->SetBinContent(i, -pulse_shape->GetBinContent(i));}
 
494
 
 
495
    for (int i = 0; i < nADC; i++) {
 
496
      pulse_shape->Fill(_rand->Landau(arrival_time, _pulse_shape_landau_width));}
 
497
 
 
498
    for (int i = 1; i <= _acquisition_window; i++) {
 
499
      pulse_shape->SetBinContent(i, baseline-pulse_shape->GetBinContent(i));}
 
500
 
 
501
    double pedestal_baseline = pulse_shape->Integral(1, _signal_integration_window)
 
502
                             /static_cast<double>(_signal_integration_window);
 
503
    double pedestal_area = pedestal_baseline*static_cast<double>(_signal_integration_window)
 
504
                         - pulse_shape->Integral(_signal_integration_window+1,
 
505
                           _signal_integration_window*2);
 
506
    double signal_area = pedestal_baseline*static_cast<double>(_signal_integration_window)
 
507
                       - pulse_shape->Integral(arrival_time-10+1,
 
508
                                               arrival_time+_signal_integration_window-10);
 
509
 
 
510
    // set pulse area
 
511
    double xAreaDigi = signal_area;
 
512
 
 
513
    // correct non physical values
 
514
    if (nADC == 0) xAreaDigi = 0;
 
515
    if (xAreaDigi < 0) xAreaDigi = 0;
 
516
 
 
517
    // set pedestal area
 
518
    double xPedestalAreaDigi = pedestal_area;
 
519
 
 
520
    // set hit arrival time
 
521
    int xArrivalTimeDigi = arrival_time;
 
522
 
 
523
    // simulate time delay in cables
 
524
    // !!! TO DO !!!
 
525
 
 
526
    // set pulse shape
 
527
    std::vector<int> xSamplesDigi;
 
528
    for (Int_t i = 1; i <= _acquisition_window; i++) {
 
529
      xSamplesDigi.push_back(pulse_shape->GetBinContent(i));}
 
530
 
 
531
    fADCdata data;
 
532
    data._orientation = xOri;
 
533
    data._charge = xAreaDigi;/*ADC*/
 
534
    data._pedestal_area = xPedestalAreaDigi;/*ADC*/
 
535
    data._time = xArrivalTimeDigi;/*ADC*/
 
536
    data._samples = xSamplesDigi;/*#*/
 
537
    emr_fadc_events_tmp[xPe][xPlane] = data;
 
538
 
 
539
    delete pulse_shape;
 
540
  }
 
541
}
 
542
 
 
543
void MapCppEMRMCDigitization::fill(MAUS::Spill *spill,
 
544
                                   int nPartEvents,
 
545
                                   EMRDBBEventVector emr_dbb_events_tmp,
 
546
                                   EMRfADCEventVector emr_fadc_events_tmp) const {
 
547
 
 
548
  int recPartEvents = spill->GetReconEventSize();
 
549
  ReconEventPArray *recEvts =  spill->GetReconEvents();
 
550
  int xSpill = spill->GetSpillNumber();
 
551
 
 
552
  // Resize the recon event to harbour all the EMR noise+decays
 
553
  if (recPartEvents == 0) { // No recEvts yet
 
554
      for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
555
        recEvts->push_back(new ReconEvent);
 
556
    }
 
557
  }
 
558
 
 
559
  if (nPartEvents == recPartEvents+2) { // Regular sized recEvts already created
 
560
    for (int iPe = 0; iPe < 2; iPe++) {
 
561
      recEvts->push_back(new ReconEvent);
 
562
    }
 
563
  }
 
564
 
 
565
  // Fill it with DBB and fADC arrays
 
566
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
567
    EMREvent *evt = new EMREvent;
 
568
    EMRPlaneHitArray plArray;
 
569
 
 
570
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
571
      int xOri  = emr_fadc_events_tmp[iPe][iPlane]._orientation;
 
572
      double xCharge = emr_fadc_events_tmp[iPe][iPlane]._charge;
 
573
      int xPos = emr_fadc_events_tmp[iPe][iPlane]._time;
 
574
      double xPedestalArea = emr_fadc_events_tmp[iPe][iPlane]._pedestal_area;
 
575
      std::vector<int> xSamples = emr_fadc_events_tmp[iPe][iPlane]._samples;
 
576
 
 
577
      EMRPlaneHit *plHit = new EMRPlaneHit;
 
578
      plHit->SetPlane(iPlane);
 
579
      plHit->SetTrigger(iPe);
 
580
      plHit->SetRun(0);
 
581
      plHit->SetOrientation(xOri);
 
582
      plHit->SetCharge(xCharge);
 
583
      plHit->SetDeltaT(xPos);
 
584
      plHit->SetSpill(xSpill);
 
585
      plHit->SetPedestalArea(xPedestalArea);
 
586
      plHit->SetSamples(xSamples);
 
587
 
 
588
      EMRBarArray barArray;
 
589
      for (int iBar = 0; iBar < _number_of_bars; iBar++) {
 
590
        int nHits = emr_dbb_events_tmp[iPe][iPlane][iBar].size();
 
591
        if (nHits) {
 
592
          EMRBar *bar = new EMRBar;
 
593
          bar->SetBar(iBar);
 
594
          bar->SetEMRBarHitArray(emr_dbb_events_tmp[iPe][iPlane][iBar]);
 
595
          barArray.push_back(bar);
 
596
        }
 
597
      }
 
598
 
 
599
      plHit->SetEMRBarArray(barArray);
 
600
      if (barArray.size() || xCharge) {
 
601
        plArray.push_back(plHit);
 
602
      }
 
603
    }
 
604
 
 
605
    if (iPe < 1) {evt->SetInitialTrigger(true);
 
606
    } else {evt->SetInitialTrigger(false);}
 
607
 
 
608
    evt->SetEMRPlaneHitArray(plArray);
 
609
    recEvts->at(iPe)->SetEMREvent(evt);
 
610
    recEvts->at(iPe)->SetPartEventNumber(iPe);
 
611
  }
 
612
 
 
613
  spill->SetReconEvents(recEvts);
 
614
}
 
615
 
 
616
EMRDBBEventVector MapCppEMRMCDigitization::get_dbb_data_tmp(int nPartEvts) const {
 
617
  EMRDBBEventVector emr_dbb_events_tmp;
 
618
  emr_dbb_events_tmp.resize(nPartEvts);
 
619
  for (int iPe = 0; iPe < nPartEvts; iPe++) {
 
620
    emr_dbb_events_tmp[iPe].resize(_number_of_planes);  // number of planes
 
621
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
622
      emr_dbb_events_tmp[iPe][iPlane].resize(_number_of_bars); // number of bars in a plane
 
623
    }
 
624
  }
 
625
  return emr_dbb_events_tmp;
 
626
}
 
627
 
 
628
EMRfADCEventVector MapCppEMRMCDigitization::get_fadc_data_tmp(int nPartEvts) const {
 
629
  EMRfADCEventVector emr_fadc_events_tmp;
 
630
  emr_fadc_events_tmp.resize(nPartEvts);
 
631
  for (int iPe = 0; iPe < nPartEvts ;iPe++) {
 
632
    emr_fadc_events_tmp[iPe].resize(_number_of_planes);
 
633
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
634
      fADCdata data;
 
635
      data._orientation = iPlane%2;
 
636
      data._charge = 0.0;
 
637
      data._pedestal_area = 0.0;
 
638
      data._time = 0;
 
639
      std::vector<int> xSamples;
 
640
      data._samples = xSamples;
 
641
      emr_fadc_events_tmp[iPe][iPlane] = data;
 
642
    }
 
643
  }
 
644
  return emr_fadc_events_tmp;
 
645
}
 
646
} // Namespace MAUS