← Back to branch summary

~durga/maus/rel709

« back to all changes in this revision

Viewing changes to src/map/MapCppEMRRecon/MapCppEMRRecon.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:
src/map/MapCppEMRRecon/MapCppEMRRecon.cc
 
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 "Utils/DAQChannelMap.hh"
 
21
#include "Interface/dataCards.hh"
 
22
#include "API/PyWrapMapBase.hh"
 
23
#include "Converter/DataConverters/CppJsonSpillConverter.hh"
 
24
#include "Converter/DataConverters/JsonCppSpillConverter.hh"
 
25
 
 
26
#include "src/map/MapCppEMRRecon/MapCppEMRRecon.hh"
 
27
 
 
28
namespace MAUS {
 
29
 
 
30
PyMODINIT_FUNC init_MapCppEMRRecon(void) {
 
31
  PyWrapMapBase<MAUS::MapCppEMRRecon>::PyWrapMapBaseModInit
 
32
                                                ("MapCppEMRRecon", "", "", "", "");
 
33
}
 
34
 
 
35
MapCppEMRRecon::MapCppEMRRecon()
 
36
    : MapBase<MAUS::Data>("MapCppEMRRecon") {
 
37
}
 
38
 
 
39
////////////////////////////////////////////////////////////////////////////
 
40
void MapCppEMRRecon::_birth(const std::string& argJsonConfigDocument) {
 
41
  _classname = "MapCppEMRRecon";
 
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
                      "MapCppEMRRecon::birth");
 
48
  }
 
49
 
 
50
  //  JsonCpp setup
 
51
  Json::Value configJSON;
 
52
 
 
53
  configJSON = JsonWrapper::StringToJson(argJsonConfigDocument);
 
54
 
 
55
  // Fetch variables
 
56
  _number_of_planes = configJSON["EMRnumberOfPlanes"].asInt();
 
57
  _number_of_bars = configJSON["EMRnumberOfBars"].asInt();
 
58
 
 
59
  _secondary_hits_bunching_distance = configJSON["EMRsecondaryHitsBunchingDistance"].asInt();
 
60
  _secondary_hits_bunching_width = configJSON["EMRsecondaryHitsBunchingWidth"].asInt();
 
61
 
 
62
  _primary_trigger_minXhits = configJSON["EMRprimaryTriggerMinXhits"].asInt();
 
63
  _primary_trigger_minYhits = configJSON["EMRprimaryTriggerMinYhits"].asInt();
 
64
  _primary_trigger_minNhits = configJSON["EMRprimaryTriggerMinNhits"].asInt();
 
65
  _secondary_trigger_minXhits = configJSON["EMRsecondaryTriggerMinXhits"].asInt();
 
66
  _secondary_trigger_minYhits = configJSON["EMRsecondaryTriggerMinYhits"].asInt();
 
67
  _secondary_trigger_minNhits = configJSON["EMRsecondaryTriggerMinNhits"].asInt();
 
68
  _secondary_trigger_minTot = configJSON["EMRsecondaryTriggerMinTot"].asInt();
 
69
 
 
70
  _max_secondary_to_primary_track_distance
 
71
        = configJSON["EMRmaxSecondaryToPrimaryTrackDistance"].asInt();
 
72
}
 
73
 
 
74
////////////////////////////////////////////////////////////////////////////
 
75
void MapCppEMRRecon::_death() {
 
76
}
 
77
 
 
78
////////////////////////////////////////////////////////////////////////////
 
79
void MapCppEMRRecon::_process(Data *data) const {
 
80
 
 
81
  // Get spill, break if there's no recon data
 
82
  Spill *spill = data->GetSpill();
 
83
 
 
84
  int nPartEvents = spill->GetReconEventSize();
 
85
  if (!nPartEvents)
 
86
      return;
 
87
 
 
88
  // Create DBB and fADC arrays with n+2 events (1 per trigger + noise + decays)
 
89
  EMRDBBEventVector emr_dbb_events_tmp = get_dbb_data_tmp(nPartEvents);
 
90
  EMRfADCEventVector emr_fadc_events_tmp = get_fadc_data_tmp(nPartEvents);
 
91
 
 
92
  // Create DBB and fADC arrays to future host n+1+n' events (1 per trigger + noise + 1 per decay)
 
93
  EMRDBBEventVector emr_dbb_events[3]; // preselected, primary, secondary
 
94
  EMRfADCEventVector emr_fadc_events;
 
95
  EMRTrackEventVector emr_track_events;
 
96
 
 
97
  // Fill temporary array with preselected events
 
98
  process_preselected_events(spill, emr_dbb_events_tmp, emr_fadc_events_tmp);
 
99
 
 
100
  // Bunch secondary hits and place them in new reconEvents
 
101
  process_secondary_events(emr_dbb_events_tmp, emr_fadc_events_tmp,
 
102
                           emr_dbb_events, emr_fadc_events, emr_track_events);
 
103
 
 
104
  // Fill primary and secondary events array with the most significant bar of each plane
 
105
  tot_cleaning(nPartEvents, emr_dbb_events, emr_fadc_events, emr_track_events);
 
106
 
 
107
  // Reconstruct the coordinates of each Hit
 
108
  coordinates_reconstruction(nPartEvents, emr_dbb_events, emr_fadc_events, emr_track_events);
 
109
 
 
110
  // Match the primary tracks with their decay
 
111
  track_matching(nPartEvents, emr_dbb_events, emr_fadc_events, emr_track_events);
 
112
 
 
113
  // Fill the Recon event array with Spill information (1 per trigger + noise + decays)
 
114
  fill(spill, nPartEvents, emr_dbb_events, emr_fadc_events, emr_track_events);
 
115
}
 
116
 
 
117
////////////////////////////////////////////////////////////////////////////
 
118
void MapCppEMRRecon::process_preselected_events(MAUS::Spill *spill,
 
119
                                                EMRDBBEventVector& emr_dbb_events_tmp,
 
120
                                                EMRfADCEventVector& emr_fadc_events_tmp) const {
 
121
 
 
122
  int nPartEvents = spill->GetReconEvents()->size();
 
123
 
 
124
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
125
 
 
126
    EMREvent *evt = spill->GetReconEvents()->at(iPe)->GetEMREvent();
 
127
 
 
128
    int nPlHits = evt->GetEMRPlaneHitArray().size();
 
129
 
 
130
    if ( !nPlHits ) continue;
 
131
 
 
132
    // Total amount of bars hit in one trigger
 
133
    int nBarsTotal = 0;
 
134
    for (int iPlane = 0; iPlane < nPlHits; iPlane++) {
 
135
      EMRPlaneHit *plHit = evt->GetEMRPlaneHitArray().at(iPlane);
 
136
      nBarsTotal += plHit->GetEMRBarArray().size();
 
137
    }
 
138
 
 
139
    // Fill collection of pre-selected events
 
140
    for (int iPlane = 0; iPlane < nPlHits; iPlane++) {
 
141
 
 
142
      EMRPlaneHit *plHit = evt->GetEMRPlaneHitArray().at(iPlane);
 
143
      int xPlane = plHit->GetPlane();
 
144
      int nBars = plHit->GetEMRBarArray().size();
 
145
 
 
146
      for (int iBar = 0; iBar < nBars; iBar++) {
 
147
 
 
148
        EMRBar *bar = plHit->GetEMRBarArray().at(iBar);
 
149
        int xBar = bar->GetBar();
 
150
 
 
151
        // Skip calibration channel
 
152
        if (xBar == 0) continue;
 
153
 
 
154
        int nBarHits = bar->GetEMRBarHitArray().size();
 
155
 
 
156
        for (int iBarHit = 0; iBarHit < nBarHits; iBarHit++) {
 
157
 
 
158
          EMRBarHit barHit = bar->GetEMRBarHitArray().at(iBarHit);
 
159
          int xTot  = barHit.GetTot();
 
160
          int xDeltaT = barHit.GetDeltaT();
 
161
          int xHitTime = barHit.GetHitTime();
 
162
          double x = 0.0;
 
163
          double y = 0.0;
 
164
          double z = 0.0;
 
165
 
 
166
          EMRBarHit bHit;
 
167
          bHit.SetTot(xTot);
 
168
          bHit.SetDeltaT(xDeltaT);
 
169
          bHit.SetHitTime(xHitTime);
 
170
          bHit.SetX(x);
 
171
          bHit.SetY(y);
 
172
          bHit.SetZ(z);
 
173
 
 
174
          emr_dbb_events_tmp[iPe][xPlane][xBar].push_back(bHit);
 
175
        }
 
176
      }
 
177
 
 
178
      int xOri = plHit->GetOrientation();
 
179
      int xArrivalTime = plHit->GetDeltaT();
 
180
      double xCharge = plHit->GetCharge();
 
181
      double xPedestalArea = plHit->GetPedestalArea();
 
182
      std::vector<int> xSamples = plHit->GetSamples();
 
183
 
 
184
      fADCdata data;
 
185
      data._orientation = xOri;
 
186
      data._charge = xCharge;
 
187
      data._pedestal_area = xPedestalArea;
 
188
      data._time = xArrivalTime;
 
189
      data._samples = xSamples;
 
190
      emr_fadc_events_tmp[iPe][xPlane] = data;
 
191
    }
 
192
  }
 
193
}
 
194
 
 
195
////////////////////////////////////////////////////////////////////////////
 
196
void MapCppEMRRecon::process_secondary_events(EMRDBBEventVector emr_dbb_events_tmp,
 
197
                                             EMRfADCEventVector emr_fadc_events_tmp,
 
198
                                             EMRDBBEventVector *emr_dbb_events,
 
199
                                             EMRfADCEventVector& emr_fadc_events,
 
200
                                             EMRTrackEventVector& emr_track_events) const {
 
201
 
 
202
  int nPartEvents = emr_fadc_events_tmp.size();
 
203
 
 
204
  // The last trigger corresponds to the secondary events
 
205
  int secondaryEventId = nPartEvents-1;
 
206
 
 
207
  // Save hit time of each secondary trigger in a vector
 
208
  std::vector<int> hitTimeVector;
 
209
  for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
210
    for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
211
      int nHits = emr_dbb_events_tmp[secondaryEventId][iPlane][iBar].size();
 
212
 
 
213
      for (int iBarHit = 0; iBarHit < nHits; iBarHit++) {
 
214
        EMRBarHit bHit = emr_dbb_events_tmp[secondaryEventId][iPlane][iBar][iBarHit];
 
215
        int xTot = bHit.GetTot(); // DBB counts
 
216
        int hitTime = bHit.GetHitTime(); // DBB counts
 
217
        if (xTot > _secondary_trigger_minTot) hitTimeVector.push_back(hitTime);
 
218
      }
 
219
    }
 
220
  }
 
221
 
 
222
  // Sort the vector in ascending order
 
223
  sort(hitTimeVector.begin(), hitTimeVector.end());
 
224
  int hitTimeVectorSize = hitTimeVector.size();
 
225
 
 
226
  // Find groups of hits defined by _secondary_hits_bunching_width
 
227
  std::vector<int> hitTimeGroup;
 
228
  int nHitsInGroup = 0;
 
229
 
 
230
  for (int iHit = 1; iHit < hitTimeVectorSize+1; iHit++) {
 
231
    int deltat;
 
232
    if (iHit == hitTimeVectorSize) deltat = _secondary_hits_bunching_distance+1;
 
233
    else
 
234
      deltat = hitTimeVector.at(iHit) - hitTimeVector.at(iHit-1);
 
235
 
 
236
    if (nHitsInGroup >= _secondary_trigger_minNhits &&
 
237
        deltat > _secondary_hits_bunching_distance ) {
 
238
      hitTimeGroup.push_back(hitTimeVector.at(iHit-1));
 
239
      nHitsInGroup = 0;
 
240
    } else if (deltat <= _secondary_hits_bunching_width) {
 
241
      nHitsInGroup++;
 
242
    }
 
243
  }
 
244
 
 
245
  // Number of the secondary trigger
 
246
  int nSeconPartEvents = hitTimeGroup.size();
 
247
 
 
248
  // Resize the event arrays to accomodate one extra event per secondary track (n')
 
249
  for (int i = 0; i < 3; i++)
 
250
    emr_dbb_events[i] = get_dbb_data_tmp(nPartEvents+nSeconPartEvents);
 
251
  emr_fadc_events = get_fadc_data_tmp(nPartEvents+nSeconPartEvents);
 
252
  emr_track_events = get_track_data_tmp(nPartEvents+nSeconPartEvents);
 
253
 
 
254
  // Fill collection of events (pre-selected)
 
255
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
256
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
257
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
258
        int nHits = emr_dbb_events_tmp[iPe][iPlane][iBar].size();
 
259
 
 
260
        for (int iBarHit = 0; iBarHit < nHits; iBarHit++) {
 
261
          EMRBarHit bHit = emr_dbb_events_tmp[iPe][iPlane][iBar][iBarHit];
 
262
          emr_dbb_events[0][iPe][iPlane][iBar].push_back(bHit);
 
263
        }
 
264
      }
 
265
 
 
266
      if (iPe < nPartEvents-2) { // All triggers except the noise and secondary
 
267
        fADCdata data;
 
268
        data._orientation = emr_fadc_events_tmp[iPe][iPlane]._orientation;
 
269
        data._charge = emr_fadc_events_tmp[iPe][iPlane]._charge;
 
270
        data._pedestal_area = emr_fadc_events_tmp[iPe][iPlane]._pedestal_area;
 
271
        data._time = emr_fadc_events_tmp[iPe][iPlane]._time;
 
272
        data._samples = emr_fadc_events_tmp[iPe][iPlane]._samples;
 
273
        emr_fadc_events[iPe][iPlane] = data;
 
274
      }
 
275
    }
 
276
  }
 
277
 
 
278
  // Associate each bunch to a secondary trigger
 
279
  for (int iHitTimeGroup = 0; iHitTimeGroup < nSeconPartEvents; iHitTimeGroup++) {
 
280
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
281
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
282
        int nHits = emr_dbb_events_tmp[secondaryEventId][iPlane][iBar].size();
 
283
 
 
284
        for (int iBarHit = 0; iBarHit < nHits; iBarHit++) {
 
285
          EMRBarHit bHit = emr_dbb_events_tmp[secondaryEventId][iPlane][iBar][iBarHit];
 
286
          int hitTime = bHit.GetHitTime(); // DBB counts
 
287
          int deltat = hitTimeGroup.at(iHitTimeGroup) - hitTime;
 
288
 
 
289
          if (deltat >= 0 && deltat < _secondary_hits_bunching_width) {
 
290
            emr_dbb_events[0][secondaryEventId+1+iHitTimeGroup]
 
291
                          [iPlane][iBar].push_back(bHit);
 
292
          }
 
293
        }
 
294
      }
 
295
    }
 
296
  }
 
297
}
 
298
 
 
299
////////////////////////////////////////////////////////////////////////////
 
300
void MapCppEMRRecon::tot_cleaning(int nPartEvents,
 
301
                                  EMRDBBEventVector *emr_dbb_events,
 
302
                                  EMRfADCEventVector& emr_fadc_events,
 
303
                                  EMRTrackEventVector& emr_track_events) const {
 
304
 
 
305
  int nTotalPartEvents = emr_fadc_events.size();
 
306
 
 
307
  for (int iPe = 0; iPe < nTotalPartEvents; iPe++) {
 
308
    // Skip noise and secondary triggers
 
309
    if (iPe == nPartEvents-2 || iPe == nPartEvents-1) continue;
 
310
 
 
311
    // Count number of X and Y planes hit
 
312
    int xPlaneHits = 0;
 
313
    int yPlaneHits = 0;
 
314
 
 
315
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
316
      int nHits = 0;
 
317
 
 
318
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
319
        nHits = nHits + emr_dbb_events[0][iPe][iPlane][iBar].size();
 
320
      }
 
321
      if (nHits && iPlane % 2 == 0) xPlaneHits++;
 
322
      if (nHits && iPlane % 2 == 1) yPlaneHits++;
 
323
    }
 
324
 
 
325
    // Reject triggers with too small amount of hits
 
326
    if (iPe < nPartEvents-2 &&
 
327
        (xPlaneHits < _primary_trigger_minXhits ||
 
328
         yPlaneHits < _primary_trigger_minYhits ||
 
329
         (xPlaneHits + yPlaneHits) < _primary_trigger_minNhits)) continue;
 
330
    if (iPe >= nPartEvents &&
 
331
        (xPlaneHits < _secondary_trigger_minXhits ||
 
332
         yPlaneHits < _secondary_trigger_minYhits ||
 
333
         (xPlaneHits + yPlaneHits) < _secondary_trigger_minNhits)) continue;
 
334
 
 
335
    emr_track_events[iPe]._has_primary = true;
 
336
 
 
337
    // Select the channel with the highest ToT in each plane
 
338
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
339
      int xTotMax  = 0;
 
340
      int xBarMax = 0;
 
341
      int xHitMax = 0;
 
342
 
 
343
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
344
        int nHits = emr_dbb_events[0][iPe][iPlane][iBar].size();
 
345
 
 
346
        for (int iBarHit = 0; iBarHit < nHits; iBarHit++) {
 
347
          EMRBarHit bHit = emr_dbb_events[0][iPe][iPlane][iBar][iBarHit];
 
348
          int xTot = bHit.GetTot();
 
349
 
 
350
          if (xTot > xTotMax) {
 
351
            xTotMax = xTot;
 
352
            xHitMax = iBarHit;
 
353
            xBarMax = iBar;
 
354
          }
 
355
        }
 
356
      }
 
357
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
358
        int nHits = emr_dbb_events[0][iPe][iPlane][iBar].size();
 
359
 
 
360
        // Fill primary events array
 
361
        for (int iBarHit = 0; iBarHit < nHits; iBarHit++) {
 
362
          EMRBarHit bHit = emr_dbb_events[0][iPe][iPlane][iBar][iBarHit];
 
363
 
 
364
          if (iBar == xBarMax && iBarHit == xHitMax &&
 
365
              bHit.GetTot()>_secondary_trigger_minTot) {
 
366
            emr_dbb_events[1][iPe][iPlane][iBar].push_back(bHit);
 
367
          }
 
368
        }
 
369
      }
 
370
    }
 
371
  }
 
372
}
 
373
 
 
374
////////////////////////////////////////////////////////////////////////////
 
375
void MapCppEMRRecon::coordinates_reconstruction(int nPartEvents,
 
376
                                                EMRDBBEventVector *emr_dbb_events,
 
377
                                                EMRfADCEventVector& emr_fadc_events,
 
378
                                                EMRTrackEventVector& emr_track_events) const {
 
379
 
 
380
  int nTotalPartEvents = emr_fadc_events.size();
 
381
 
 
382
  for (int iPe = 0; iPe < nTotalPartEvents; iPe++) {
 
383
 
 
384
    // Skip noise and secondary triggers
 
385
    if (iPe == nPartEvents-1 || iPe == nPartEvents-2) continue;
 
386
 
 
387
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
388
 
 
389
      EMRBarHit Hit0;
 
390
      bool Hit0Found = false;
 
391
      bool Hit1Found = false;
 
392
      bool Hit2Found = false;
 
393
 
 
394
      int barid = -1;
 
395
 
 
396
      int x0(-1), x1(-1), x2(-1), y1(-1), y2(-1);
 
397
      double y0(-1.0), a(-1.0), b(-1.0);
 
398
 
 
399
      // Find Primary hit
 
400
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
401
        if (Hit0Found) break;
 
402
        if (emr_dbb_events[1][iPe][iPlane][iBar].size()) {
 
403
          Hit0 = emr_dbb_events[1][iPe][iPlane][iBar][0];
 
404
          Hit0Found = true;
 
405
          x0 = iPlane;
 
406
          barid = iBar;
 
407
        }
 
408
      }
 
409
 
 
410
      if (!Hit0Found) continue;
 
411
 
 
412
      // Look backwards for hits
 
413
      for (int aPlane = iPlane-1; aPlane >= 0; aPlane = aPlane-2) {
 
414
        if (Hit1Found) break;
 
415
        for (int aBar = 1; aBar < _number_of_bars; aBar++) { // Skip test channel
 
416
          if (emr_dbb_events[1][iPe][aPlane][aBar].size()) {
 
417
            if (!Hit1Found) {
 
418
              Hit1Found = true;
 
419
              x1 = aPlane;
 
420
              y1 = aBar;
 
421
            }
 
422
          }
 
423
        }
 
424
      }
 
425
 
 
426
      // Look forward for hits
 
427
      for (int bPlane = iPlane+1; bPlane < _number_of_planes; bPlane = bPlane+2) {
 
428
        if (Hit1Found && Hit2Found) break;
 
429
        for (int bBar = 1; bBar < _number_of_bars; bBar++) { // Skip test channel
 
430
          if (emr_dbb_events[1][iPe][bPlane][bBar].size()) {
 
431
            if (Hit2Found && !Hit1Found) {
 
432
              Hit1Found = true;
 
433
              x1 = bPlane;
 
434
              y1 = bBar;
 
435
            }
 
436
            if (!Hit2Found) {
 
437
              Hit2Found = true;
 
438
              x2 = bPlane;
 
439
              y2 = bBar;
 
440
            }
 
441
          }
 
442
        }
 
443
      }
 
444
 
 
445
      // Look backwards for the second hit if nothing found in the forward direction
 
446
      for (int aPlane = iPlane-1; aPlane >= 0; aPlane = aPlane-2) {
 
447
        if (Hit1Found && Hit2Found) break;
 
448
        if (aPlane == x1) continue;
 
449
        for (int aBar = 1; aBar < _number_of_bars; aBar++) { // Skip test channel
 
450
          if (emr_dbb_events[1][iPe][aPlane][aBar].size()) {
 
451
            if (Hit1Found && !Hit2Found) {
 
452
              Hit2Found = true;
 
453
              x2 = aPlane;
 
454
              y2 = aBar;
 
455
            }
 
456
          }
 
457
        }
 
458
      }
 
459
 
 
460
      // Calculate the coordinates of the hit in each plane
 
461
      if (Hit1Found && Hit2Found) {
 
462
        a = (static_cast<double>(y2)-static_cast<double>(y1))
 
463
            / (static_cast<double>(x2)-static_cast<double>(x1));
 
464
        b = static_cast<double>(y1) - a*static_cast<double>(x1);
 
465
        y0 = a*static_cast<double>(x0) + b;
 
466
      }
 
467
 
 
468
      if (Hit1Found && !Hit2Found) {
 
469
        y0 = y1;
 
470
      }
 
471
 
 
472
      if (!Hit1Found && Hit2Found) {
 
473
        y0 = y2;
 
474
      }
 
475
 
 
476
      if (iPlane % 2 == 0) {
 
477
        emr_dbb_events[1][iPe][x0][barid][0].SetX(barid);
 
478
        emr_dbb_events[1][iPe][x0][barid][0].SetY(y0);
 
479
      }
 
480
      if (iPlane % 2 == 1) {
 
481
        emr_dbb_events[1][iPe][x0][barid][0].SetX(y0);
 
482
        emr_dbb_events[1][iPe][x0][barid][0].SetY(barid);
 
483
      }
 
484
      emr_dbb_events[1][iPe][x0][barid][0].SetZ(x0);
 
485
    }
 
486
  }
 
487
}
 
488
 
 
489
////////////////////////////////////////////////////////////////////////////
 
490
void MapCppEMRRecon::track_matching(int nPartEvents,
 
491
                                    EMRDBBEventVector *emr_dbb_events,
 
492
                                    EMRfADCEventVector& emr_fadc_events,
 
493
                                    EMRTrackEventVector& emr_track_events) const {
 
494
 
 
495
  int nTotalPartEvents = emr_fadc_events.size();
 
496
 
 
497
  for (int iPe = 0; iPe < nPartEvents; iPe++) {
 
498
    // Range calculation of primary tracks
 
499
    double x1(-1.0), y1(-1.0), z1(-1.0);
 
500
    double x2(-1.0), y2(-1.0), z2(-1.0);
 
501
    double primEventRange = 0.0;
 
502
    bool primHitsFound = false;
 
503
 
 
504
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
505
      for (int iBar = 1; iBar < _number_of_bars; iBar++) { // Skip test channel
 
506
        if (emr_dbb_events[1][iPe][iPlane][iBar].size()) {
 
507
          EMRBarHit bHit = emr_dbb_events[1][iPe][iPlane][iBar][0];
 
508
          // Last point
 
509
          x2 = bHit.GetX();
 
510
          y2 = bHit.GetY();
 
511
          z2 = bHit.GetZ();
 
512
          if (primHitsFound) {
 
513
            primEventRange = primEventRange + sqrt(pow(x1-x2, 2)+pow(y1-y2, 2)+pow(z1-z2, 2));
 
514
          }
 
515
          // Previous point
 
516
          x1 = bHit.GetX();
 
517
          y1 = bHit.GetY();
 
518
          z1 = bHit.GetZ();
 
519
          primHitsFound = true;
 
520
        }
 
521
      }
 
522
    }
 
523
 
 
524
    // Range calculation of secondary tracks
 
525
    double minDist = 999.9;
 
526
    int seconEventId = -1;
 
527
    double seconEventRange = -1.0;
 
528
 
 
529
    for (int iiPe = nPartEvents; iiPe < nTotalPartEvents; iiPe++) {
 
530
      double x3(-1.0), y3(-1.0), z3(-1.0);
 
531
      double x4(-1.0), y4(-1.0), z4(-1.0);
 
532
      double x5(-1.0), y5(-1.0), z5(-1.0);
 
533
      double seconEventRangeTmp = 0.0;
 
534
      bool seconHitsFound = false;
 
535
 
 
536
      for (int iiPlane = 0; iiPlane < _number_of_planes; iiPlane++) {
 
537
        for (int iiBar = 1; iiBar < _number_of_bars; iiBar++) { // Skip test channel
 
538
          if (emr_dbb_events[1][iiPe][iiPlane][iiBar].size()) {
 
539
            EMRBarHit bHit = emr_dbb_events[1][iiPe][iiPlane][iiBar][0];
 
540
 
 
541
            // Last point
 
542
            x5 = bHit.GetX();
 
543
            y5 = bHit.GetY();
 
544
            z5 = bHit.GetZ();
 
545
 
 
546
            if (!seconHitsFound) {
 
547
              // First point
 
548
              x3 = x5;
 
549
              y3 = y5;
 
550
              z3 = z5;
 
551
            }
 
552
            if (seconHitsFound) {
 
553
              seconEventRangeTmp = seconEventRangeTmp
 
554
                                 + sqrt(pow(x5-x4, 2)+pow(y5-y4, 2)+pow(z5-z4, 2));
 
555
            }
 
556
            // Previous point
 
557
            x4 = bHit.GetX();
 
558
            y4 = bHit.GetY();
 
559
            z4 = bHit.GetZ();
 
560
            seconHitsFound = true;
 
561
          }
 
562
        }
 
563
      }
 
564
 
 
565
      // Look for a track close enough to the primary to match it
 
566
      if (seconHitsFound) {
 
567
        double dist23 = sqrt(pow(x2-x3, 2)+pow(y2-y3, 2)+pow(z2-z3, 2));
 
568
        double dist25 = sqrt(pow(x2-x5, 2)+pow(y2-y5, 2)+pow(z2-z5, 2));
 
569
        double dist = std::min(dist23, dist25);
 
570
 
 
571
        if (dist < minDist && dist < _max_secondary_to_primary_track_distance) {
 
572
          minDist = dist;
 
573
          seconEventId = iiPe;
 
574
          seconEventRange = seconEventRangeTmp;
 
575
        }
 
576
      }
 
577
    }
 
578
 
 
579
//    std::cerr << "minDist: " << minDist << " seconEventId: " << seconEventId << std::endl;
 
580
//    std::cerr << "Rp=" << primEventRange << " Rs=" << seconEventRange << "   " << std::endl;
 
581
 
 
582
    // If a secondary is found, fill the corresponding secondary array
 
583
    if (seconEventId != -1) {
 
584
      for (int iiPlane = 0; iiPlane < _number_of_planes; iiPlane++) {
 
585
        for (int iiBar = 1; iiBar < _number_of_bars; iiBar++) { // Skip test channel
 
586
          if (emr_dbb_events[1][seconEventId][iiPlane][iiBar].size()) {
 
587
            EMRBarHit bHit = emr_dbb_events[1][seconEventId][iiPlane][iiBar][0];
 
588
            emr_dbb_events[2][iPe][iiPlane][iiBar].push_back(bHit);
 
589
          }
 
590
        }
 
591
      }
 
592
      emr_track_events[iPe]._range_secondary = seconEventRange;
 
593
      emr_track_events[iPe]._has_secondary = true;
 
594
      emr_track_events[iPe]._secondary_to_primary_track_distance = minDist;
 
595
    }
 
596
    emr_track_events[iPe]._range_primary = primEventRange;
 
597
  }
 
598
}
 
599
 
 
600
////////////////////////////////////////////////////////////////////////////
 
601
void MapCppEMRRecon::fill(Spill *spill,
 
602
                          int nPartEvents,
 
603
                          EMRDBBEventVector *emr_dbb_events,
 
604
                          EMRfADCEventVector& emr_fadc_events,
 
605
                          EMRTrackEventVector& emr_track_events) const {
 
606
 
 
607
  int nTotalPartEvents = emr_fadc_events.size();
 
608
  int recPartEvents = spill->GetReconEventSize();
 
609
  int xRun = spill->GetRunNumber();
 
610
  int xSpill = spill->GetSpillNumber();
 
611
 
 
612
  ReconEventPArray *recEvts =  spill->GetReconEvents();
 
613
 
 
614
  if (recPartEvents < nTotalPartEvents) {
 
615
    for (int iPe = recPartEvents; iPe < nTotalPartEvents; iPe++)
 
616
      recEvts->push_back(new ReconEvent);
 
617
  }
 
618
 
 
619
  for (int iPe = 0; iPe < nTotalPartEvents; iPe++) {
 
620
    EMREvent *evt = new EMREvent;
 
621
    EMRPlaneHitArray plArray;
 
622
 
 
623
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
624
      int xOri  = emr_fadc_events[iPe][iPlane]._orientation;
 
625
      double xCharge = emr_fadc_events[iPe][iPlane]._charge;
 
626
      int xArrivalTime = emr_fadc_events[iPe][iPlane]._time;
 
627
      double xPedestalArea = emr_fadc_events[iPe][iPlane]._pedestal_area;
 
628
      std::vector<int> xSamples = emr_fadc_events[iPe][iPlane]._samples;
 
629
 
 
630
      EMRPlaneHit *plHit = new EMRPlaneHit;
 
631
      plHit->SetPlane(iPlane);
 
632
      plHit->SetTrigger(iPe);
 
633
      plHit->SetRun(xRun);
 
634
      plHit->SetOrientation(xOri);
 
635
      plHit->SetCharge(xCharge);
 
636
      plHit->SetDeltaT(xArrivalTime);
 
637
      plHit->SetSpill(xSpill);
 
638
      plHit->SetPedestalArea(xPedestalArea);
 
639
      plHit->SetSamples(xSamples);
 
640
 
 
641
      EMRBarArray barArray;
 
642
      EMRBarArray barArrayPrimary;
 
643
      EMRBarArray barArraySecondary;
 
644
      for (int i = 0; i < 3; i++) {
 
645
        for (int iBar = 1; iBar < _number_of_bars; iBar++) {
 
646
          int nHits = emr_dbb_events[i][iPe][iPlane][iBar].size();
 
647
          if (nHits) {
 
648
            EMRBar *bar = new EMRBar;
 
649
            bar->SetBar(iBar);
 
650
            bar->SetEMRBarHitArray(emr_dbb_events[i][iPe][iPlane][iBar]);
 
651
            if (i == 0) barArray.push_back(bar);
 
652
            if (i == 1) barArrayPrimary.push_back(bar);
 
653
            if (i == 2) barArraySecondary.push_back(bar);
 
654
          }
 
655
        }
 
656
      }
 
657
 
 
658
      plHit->SetEMRBarArray(barArray);
 
659
      plHit->SetEMRBarArrayPrimary(barArrayPrimary);
 
660
      plHit->SetEMRBarArraySecondary(barArraySecondary);
 
661
      if (barArray.size() || barArrayPrimary.size() || barArraySecondary.size() || xCharge) {
 
662
        plArray.push_back(plHit);
 
663
      }
 
664
    }
 
665
    evt->SetEMRPlaneHitArray(plArray);
 
666
    evt->SetRangePrimary(emr_track_events[iPe]._range_primary);
 
667
    evt->SetRangeSecondary(emr_track_events[iPe]._range_secondary);
 
668
    evt->SetSecondaryToPrimaryTrackDistance
 
669
         (emr_track_events[iPe]._secondary_to_primary_track_distance);
 
670
    evt->SetHasPrimary(emr_track_events[iPe]._has_primary);
 
671
    evt->SetHasSecondary(emr_track_events[iPe]._has_secondary);
 
672
 
 
673
    if (iPe < nPartEvents-2) evt->SetInitialTrigger(true);
 
674
    else
 
675
      evt->SetInitialTrigger(false);
 
676
 
 
677
    // std::cerr << "************************************************" << std::endl;
 
678
    // std::cerr << "range_primary = " << emr_track_events[iPe]._range_primary << std::endl;
 
679
    // std::cerr << "range_secondary = " << emr_track_events[iPe]._range_secondary << std::endl;
 
680
    // std::cerr << "secondary_to_primary_track_distance = "
 
681
    //           << emr_track_events[iPe]._secondary_to_primary_track_distance << std::endl;
 
682
    // std::cerr << "has_primary = " << emr_track_events[iPe]._has_primary << std::endl;
 
683
    // std::cerr << "has_secondary = " << emr_track_events[iPe]._has_secondary << std::endl;
 
684
    // std::cerr << "================================================" << std::endl;
 
685
 
 
686
    recEvts->at(iPe)->SetEMREvent(evt);
 
687
    recEvts->at(iPe)->SetPartEventNumber(iPe);
 
688
  }
 
689
  spill->SetReconEvents(recEvts);
 
690
}
 
691
 
 
692
EMRDBBEventVector MapCppEMRRecon::get_dbb_data_tmp(int nPartEvts) const {
 
693
  EMRDBBEventVector emr_dbb_events_tmp;
 
694
  emr_dbb_events_tmp.resize(nPartEvts);
 
695
  for (int iPe = 0; iPe < nPartEvts; iPe++) {
 
696
    emr_dbb_events_tmp[iPe].resize(_number_of_planes);  // number of planes
 
697
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
698
      emr_dbb_events_tmp[iPe][iPlane].resize(_number_of_bars); // number of bars in a plane
 
699
    }
 
700
  }
 
701
  return emr_dbb_events_tmp;
 
702
}
 
703
 
 
704
EMRfADCEventVector MapCppEMRRecon::get_fadc_data_tmp(int nPartEvts) const {
 
705
  EMRfADCEventVector emr_fadc_events_tmp;
 
706
  emr_fadc_events_tmp.resize(nPartEvts);
 
707
  for (int iPe = 0; iPe < nPartEvts ;iPe++) {
 
708
    emr_fadc_events_tmp[iPe].resize(_number_of_planes);
 
709
    for (int iPlane = 0; iPlane < _number_of_planes; iPlane++) {
 
710
      fADCdata data;
 
711
      data._orientation = iPlane%2;
 
712
      data._charge = 0.0;
 
713
      data._pedestal_area = 0.0;
 
714
      data._time = 0;
 
715
      std::vector<int> xSamples;
 
716
      data._samples = xSamples;
 
717
      emr_fadc_events_tmp[iPe][iPlane] = data;
 
718
    }
 
719
  }
 
720
  return emr_fadc_events_tmp;
 
721
}
 
722
 
 
723
EMRTrackEventVector MapCppEMRRecon::get_track_data_tmp(int nPartEvts) const {
 
724
  EMRTrackEventVector emr_track_events_tmp;
 
725
  emr_track_events_tmp.resize(nPartEvts);
 
726
  for (int iPe = 0; iPe < nPartEvts ;iPe++) {
 
727
    TrackData data;
 
728
    data._range_primary = 0.0;
 
729
    data._range_secondary = 0.0;
 
730
    data._secondary_to_primary_track_distance = 0.0;
 
731
    data._has_primary = false;
 
732
    data._has_secondary = false;
 
733
    emr_track_events_tmp[iPe] = data;
 
734
  }
 
735
  return emr_track_events_tmp;
 
736
}
 
737
}