← Back to branch summary

~maus-release/maus/release

« back to all changes in this revision

Viewing changes to src/reduce/ReduceCppEVExport/EVExporter.cc

  • Committer: Paolo Franchini
  • Date: 2018-06-24 14:27:26 UTC
  • mfrom: (659.2.80 release-candidate)
  • Revision ID: p.franchini@warwick.ac.uk-20180624142726-nbxxyjyer146dr1t
Tags: MAUS-v3.2.0
MAUS-v3.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/reduce/ReduceCppEVExport/EVExporter.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 "src/reduce/ReduceCppEVExport/EVExporter.hh"
 
19
 
 
20
#include <vector>
 
21
 
 
22
#include "DataStructure/SciFiTrack.hh"
 
23
#include "DataStructure/SciFiTrackPoint.hh"
 
24
 
 
25
namespace EventViewer {
 
26
 
 
27
EVExporter::EVExporter(MAUS::Spill *sp) {
 
28
  spill = sp;
 
29
}
 
30
 
 
31
EVExporter::~EVExporter() {
 
32
  // Do nothing
 
33
}
 
34
 
 
35
void EVExporter::ReadAllEvents(int eventSelection, bool verbose) {
 
36
  // - iterate through events
 
37
  for (size_t i = 0; i < spill->GetReconEvents()->size(); ++i) {
 
38
 
 
39
    detectorsHit = 0;
 
40
    evEvent.Reset();
 
41
    evEvent.runNumber = spill->GetRunNumber();
 
42
    evEvent.spillNumber = spill->GetSpillNumber();
 
43
    evEvent.eventNumber = i;
 
44
 
 
45
    // tof_event = new MAUS::TOFEvent; // will this help?
 
46
    tof_event = (*spill->GetReconEvents())[i]->GetTOFEvent();
 
47
    if (tof_event != NULL)
 
48
      ProcessTOFEvent();
 
49
 
 
50
    // scifi_event = new MAUS::SciFiEvent; // will this help?
 
51
    scifi_event = (*spill->GetReconEvents())[i]->GetSciFiEvent();
 
52
    if (scifi_event != NULL)
 
53
      ProcessScifiEvent();
 
54
 
 
55
    // - export one EVEvent to HepRep
 
56
    // add selection of visualization output format at some point,
 
57
    // maybe as an argument to method (only HepRep supported for now)
 
58
    if (eventSelection == 0) {
 
59
      if (verbose) {
 
60
        std::cout << "EVExporter: processing spill " << evEvent.spillNumber
 
61
                  << "   event: " << evEvent.eventNumber << std::endl;
 
62
      }
 
63
      EVHepRepExporter exp(evEvent);
 
64
      exp.Export(0);
 
65
    } else if (eventSelection == 1 && detectorsHit == 0) {
 
66
      if (verbose) {
 
67
        std::cout << "EVExporter: processing spill " << evEvent.spillNumber
 
68
                  << "   event: " << evEvent.eventNumber << std::endl;
 
69
      }
 
70
      EVHepRepExporter exp(evEvent);
 
71
      exp.Export(0);
 
72
    } else if (detectorsHit == eventSelection) {
 
73
      if (verbose) {
 
74
        std::cout << "EVExporter: processing spill " << evEvent.spillNumber
 
75
                  << "   event: " << evEvent.eventNumber << std::endl;
 
76
      }
 
77
      EVHepRepExporter exp(evEvent);
 
78
      exp.Export(0);
 
79
    }
 
80
  }
 
81
}
 
82
 
 
83
int EVExporter::ReadOneEvent(int eventNumber, std::string opt, int eventSelection, int display) {
 
84
  detectorsHit = 0;
 
85
  evEvent.Reset();
 
86
 
 
87
  evEvent.runNumber = spill->GetRunNumber();
 
88
  evEvent.spillNumber = spill->GetSpillNumber();
 
89
  evEvent.eventNumber = eventNumber;
 
90
 
 
91
  if ((*spill->GetReconEvents())[evEvent.eventNumber] == NULL) {
 
92
      std::cerr << "Event " << evEvent.eventNumber << " not available in current spill!" << "\n";
 
93
      return 0;
 
94
  }
 
95
 
 
96
  tof_event = new MAUS::TOFEvent; // will this help?
 
97
  tof_event = (*spill->GetReconEvents())[eventNumber]->GetTOFEvent();
 
98
  if (tof_event != NULL)
 
99
    ProcessTOFEvent();
 
100
 
 
101
  scifi_event = new MAUS::SciFiEvent; // will this help?
 
102
  scifi_event = (*spill->GetReconEvents())[eventNumber]->GetSciFiEvent();
 
103
  if (scifi_event != NULL)
 
104
    ProcessScifiEvent();
 
105
 
 
106
 
 
107
  // - export one EVEvent to HepRep
 
108
  // add selection of visualization output format at some point,
 
109
  // maybe as an argument to method (only HepRep supported for now)
 
110
  if (opt == "exp" && eventSelection == 0) {
 
111
    std::cout << "EVExporter: processing spill " << evEvent.spillNumber
 
112
              << "   event: " << evEvent.eventNumber << std::endl; // add vebosity option?
 
113
    EVHepRepExporter exp(evEvent);
 
114
    exp.Export(display);
 
115
  } else if (opt == "exp" && eventSelection == 1 && detectorsHit == 0) {
 
116
    std::cout << "EVExporter: processing spill " << evEvent.spillNumber
 
117
              << "   event: " << evEvent.eventNumber << std::endl; // add vebosity option?
 
118
    EVHepRepExporter exp(evEvent);
 
119
    exp.Export(display);
 
120
  } else if (opt == "exp" && detectorsHit == eventSelection) {
 
121
    std::cout << "EVExporter: processing spill " << evEvent.spillNumber
 
122
              << "   event: " << evEvent.eventNumber << std::endl; // add vebosity option?
 
123
    EVHepRepExporter exp(evEvent);
 
124
    exp.Export(display);
 
125
  } else if (opt != "exp" && opt != "no_exp") {
 
126
    std::cerr << "Error! Wrong option passed to EVExporter::ReadEvent! "
 
127
              << "Options are 'exp' to export to HepRep or 'no_exp' not to"
 
128
              << std::endl; // revise later whether this is necessary?
 
129
    return 0;
 
130
  }
 
131
 
 
132
  if (detectorsHit > 0) {
 
133
    return detectorsHit;
 
134
  } else if (detectorsHit == 0) {
 
135
    return 1; // check whether this is working correctly for non existing events
 
136
  }
 
137
}
 
138
 
 
139
 
 
140
void EVExporter::ProcessTOFEvent() {
 
141
 
 
142
  ProcessTOF0();
 
143
  ProcessTOF1();
 
144
  ProcessTOF2();
 
145
}
 
146
 
 
147
void EVExporter::ProcessTOF0() {
 
148
 
 
149
  MAUS::TOFEventSpacePoint space_points = tof_event->GetTOFEventSpacePoint();
 
150
  MAUS::TOFSpacePoint tof0_space_points;
 
151
 
 
152
  // - modify ReadOneEvent return value
 
153
  if (space_points.GetTOF0SpacePointArray().size() > 0)
 
154
    detectorsHit += 2;
 
155
 
 
156
  // 1. Loop over TOF0 space points:
 
157
  for (size_t i = 0; i < space_points.GetTOF0SpacePointArray().size(); ++i)  {
 
158
 
 
159
    tof0_space_points = space_points.GetTOF0SpacePointArray()[i];
 
160
 
 
161
    double TOF0_x = tof0_space_points.GetGlobalPosX();
 
162
    double TOF0_y = tof0_space_points.GetGlobalPosY();
 
163
    double TOF0_z = tof0_space_points.GetGlobalPosZ();
 
164
 
 
165
    MAUS::ThreeVector tof0Coordinates(TOF0_x, TOF0_y, TOF0_z);
 
166
    evEvent.tofPoints[0] = tof0Coordinates;
 
167
  }
 
168
}
 
169
 
 
170
void EVExporter::ProcessTOF1() {
 
171
  MAUS::TOFEventSpacePoint space_points = tof_event->GetTOFEventSpacePoint();
 
172
  MAUS::TOFSpacePoint tof1_space_points;
 
173
 
 
174
  // - modify ReadOneEvent return value
 
175
  if (space_points.GetTOF1SpacePointArray().size() > 0)
 
176
    detectorsHit += 4;
 
177
 
 
178
  for (size_t i = 0; i < space_points.GetTOF1SpacePointArray().size(); ++i) {
 
179
 
 
180
    tof1_space_points = space_points.GetTOF1SpacePointArray()[i];
 
181
 
 
182
    double TOF1_x = tof1_space_points.GetGlobalPosX();
 
183
    double TOF1_y = tof1_space_points.GetGlobalPosY();
 
184
    double TOF1_z = tof1_space_points.GetGlobalPosZ();
 
185
 
 
186
    MAUS::ThreeVector tof1Coordinates(TOF1_x, TOF1_y, TOF1_z);
 
187
    evEvent.tofPoints[1] = tof1Coordinates;
 
188
  }
 
189
}
 
190
 
 
191
void EVExporter::ProcessTOF2() {
 
192
  MAUS::TOFEventSpacePoint space_points = tof_event->GetTOFEventSpacePoint();
 
193
  MAUS::TOFSpacePoint tof2_space_points;
 
194
 
 
195
  // - modify ReadOneEvent return value
 
196
  if (space_points.GetTOF2SpacePointArray().size() > 0)
 
197
    detectorsHit += 32;
 
198
 
 
199
  for (size_t i = 0; i < space_points.GetTOF2SpacePointArray().size(); ++i) {
 
200
 
 
201
    tof2_space_points = space_points.GetTOF2SpacePointArray()[i];
 
202
 
 
203
    double TOF2_x = tof2_space_points.GetGlobalPosX();
 
204
    double TOF2_y = tof2_space_points.GetGlobalPosY();
 
205
    double TOF2_z = tof2_space_points.GetGlobalPosZ();
 
206
 
 
207
    MAUS::ThreeVector tof2Coordinates(TOF2_x, TOF2_y, TOF2_z);
 
208
    evEvent.tofPoints[2] = tof2Coordinates;
 
209
  }
 
210
}
 
211
 
 
212
void EVExporter::ProcessScifiEvent() {
 
213
  // Extract raw scifi cluster data
 
214
  std::vector<MAUS::SciFiCluster*> clusters = scifi_event->clusters();
 
215
  std::vector<MAUS::SciFiCluster*>::iterator clus_iter;
 
216
 
 
217
  for (clus_iter = clusters.begin(); clus_iter != clusters.end(); ++clus_iter) {
 
218
    MAUS::SciFiCluster* clus = (*clus_iter);
 
219
    int tracker = clus->get_tracker();
 
220
    int station = clus->get_station();
 
221
    int plane = clus->get_plane();
 
222
    MAUS::ThreeVector position = clus->get_position();
 
223
    double alpha = clus->get_alpha();
 
224
 
 
225
    if (tracker == 0) {
 
226
      evEvent.scifiUSTrackerClusters[plane][0].push_back(-alpha);
 
227
      evEvent.scifiUSTrackerClusters[plane][1].push_back(-position.z()); // add offset later
 
228
    } else if (tracker == 1) {
 
229
      evEvent.scifiDSTrackerClusters[plane][0].push_back(alpha);
 
230
      evEvent.scifiDSTrackerClusters[plane][1].push_back(position.z()); // add offset later
 
231
    }
 
232
  }
 
233
 
 
234
  // --- iterate through scifi tracks
 
235
  std::vector<MAUS::SciFiTrack*> tracks = scifi_event->scifitracks();
 
236
  std::vector<MAUS::SciFiTrack*>::iterator track_iter;
 
237
 
 
238
  int USPoints = 0, DSPoints = 0; // track point counters for US and DS tracker
 
239
 
 
240
  for (track_iter = tracks.begin(); track_iter != tracks.end(); ++track_iter) {
 
241
    std::vector<MAUS::SciFiTrackPoint*> track_points = (*track_iter)->scifitrackpoints();
 
242
    std::vector<MAUS::SciFiTrackPoint*>::iterator track_point_iter;
 
243
 
 
244
    // - iterate through track points
 
245
    for (track_point_iter = track_points.begin();
 
246
         track_point_iter != track_points.end(); ++track_point_iter) {
 
247
 
 
248
      MAUS::SciFiTrackPoint* point = (*track_point_iter);
 
249
      int tracker = point->tracker();
 
250
      int station = point->station();
 
251
      MAUS::ThreeVector position = point->pos();
 
252
      MAUS::ThreeVector momentum = point->mom();
 
253
      // check for nans from tracker reconstruction
 
254
      if (tracker == 0 && !(IsVectorNaN(position) || IsVectorNaN(momentum))) {
 
255
        evEvent.scifiUSTrackerPoints[station-1] = position;
 
256
        evEvent.scifiUSTrackerPointsMomenta[station-1] = momentum;
 
257
        ++USPoints;
 
258
      // check for nans from tracker reconstruction
 
259
      } else if (tracker == 1 && !(IsVectorNaN(position) || IsVectorNaN(momentum))) {
 
260
        evEvent.scifiDSTrackerPoints[station-1] = position;
 
261
        evEvent.scifiDSTrackerPointsMomenta[station-1] = momentum;
 
262
        ++DSPoints;
 
263
      }
 
264
    }
 
265
  }
 
266
 
 
267
  // - modify ReadOneEvent return value
 
268
  if (USPoints > 0)
 
269
    detectorsHit += 8;
 
270
  if (DSPoints > 0)
 
271
    detectorsHit += 16;
 
272
 
 
273
  // --- iterate through space points
 
274
  std::vector<MAUS::SciFiSpacePoint*> space_points = scifi_event->spacepoints();
 
275
  std::vector<MAUS::SciFiSpacePoint*>::iterator spoint_iter;
 
276
 
 
277
  for (spoint_iter = space_points.begin(); spoint_iter != space_points.end(); ++spoint_iter) {
 
278
    MAUS::SciFiSpacePoint* point = (*spoint_iter);
 
279
    int tracker = point->get_tracker();
 
280
 
 
281
    if (tracker == 0)
 
282
      evEvent.scifiUSTrackerSpacePoints.push_back(point->get_global_position());
 
283
    else if (tracker == 1)
 
284
      evEvent.scifiDSTrackerSpacePoints.push_back(point->get_global_position());
 
285
  }
 
286
 
 
287
  // --- iterate through straight tracks
 
288
  std::vector<MAUS::SciFiStraightPRTrack*> straight_tracks = scifi_event->straightprtracks();
 
289
  std::vector<MAUS::SciFiStraightPRTrack*>::iterator straight_track_iter;
 
290
 
 
291
  for (straight_track_iter = straight_tracks.begin();
 
292
       straight_track_iter != straight_tracks.end(); ++straight_track_iter) {
 
293
    MAUS::SciFiStraightPRTrack* strTrack = (*straight_track_iter);
 
294
 
 
295
    double x0 = strTrack->get_global_x0();
 
296
    double mx = strTrack->get_global_mx();
 
297
    double y0 = strTrack->get_global_y0();
 
298
    double my = strTrack->get_global_my();
 
299
 
 
300
    double what = strTrack->get_x0();
 
301
 
 
302
    // move these to settings instead of hardcoding them?
 
303
    double extra = 50;
 
304
    double zminUS = 13958 - extra;
 
305
    double zmaxUS = 15059 + extra;
 
306
    double zminDS = 18854 - extra;
 
307
    double zmaxDS = 19955 + extra;
 
308
 
 
309
    int tracker = strTrack->get_tracker();
 
310
    if (tracker == 0) {
 
311
      double xFirst = mx*zminUS + x0;
 
312
      double yFirst = my*zminUS + y0;
 
313
      MAUS::ThreeVector first(xFirst, yFirst, zminUS);
 
314
      evEvent.scifiUSTrackerStraightTrackPoints[0] = first;
 
315
 
 
316
      double xSecond = mx*zmaxUS + x0;
 
317
      double ySecond = my*zmaxUS + y0;
 
318
      MAUS::ThreeVector second(xSecond, ySecond, zmaxUS);
 
319
      evEvent.scifiUSTrackerStraightTrackPoints[1] = second;
 
320
    } else if (tracker == 1) {
 
321
      double xFirst = mx*zminDS + x0;
 
322
      double yFirst = my*zminDS + y0;
 
323
      MAUS::ThreeVector first(xFirst, yFirst, zminDS);
 
324
      evEvent.scifiDSTrackerStraightTrackPoints[0] = first;
 
325
 
 
326
      double xSecond = mx*zmaxDS + x0;
 
327
      double ySecond = my*zmaxDS + y0;
 
328
      MAUS::ThreeVector second(xSecond, ySecond, zmaxDS);
 
329
      evEvent.scifiDSTrackerStraightTrackPoints[1] = second;
 
330
    }
 
331
  }
 
332
 
 
333
  // --- iterate through helical tracks
 
334
  std::vector<MAUS::SciFiHelicalPRTrack*> helical_tracks = scifi_event->helicalprtracks();
 
335
  std::vector<MAUS::SciFiHelicalPRTrack*>::iterator helical_track_iter;
 
336
 
 
337
  for (helical_track_iter = helical_tracks.begin();
 
338
       helical_track_iter != helical_tracks.end(); ++helical_track_iter) {
 
339
 
 
340
      MAUS::SciFiHelicalPRTrack* heliTrack = (*helical_track_iter);
 
341
      MAUS::ThreeVector startPoint = heliTrack->get_pos0();
 
342
 
 
343
      if (heliTrack->get_tracker() == 0) {
 
344
        evEvent.scifiUSTrackerHelicalTrackParameters[0] = heliTrack->get_circle_x0();
 
345
        evEvent.scifiUSTrackerHelicalTrackParameters[1] = heliTrack->get_circle_y0();
 
346
        evEvent.scifiUSTrackerHelicalTrackParameters[2] = heliTrack->get_R();
 
347
        evEvent.scifiUSTrackerHelicalTrackParameters[3] = heliTrack->get_dsdz();
 
348
        evEvent.scifiUSTrackerHelicalTrackParameters[4] = heliTrack->get_line_sz_c();
 
349
        evEvent.scifiUSTrackerHelicalTrackParameters[5] = startPoint.x();
 
350
        evEvent.scifiUSTrackerHelicalTrackParameters[6] = startPoint.y();
 
351
        evEvent.scifiUSTrackerHelicalTrackParameters[7] = startPoint.z();
 
352
      } else if (heliTrack->get_tracker() == 1) {
 
353
        evEvent.scifiDSTrackerHelicalTrackParameters[0] = heliTrack->get_circle_x0();
 
354
        evEvent.scifiDSTrackerHelicalTrackParameters[1] = heliTrack->get_circle_y0();
 
355
        evEvent.scifiDSTrackerHelicalTrackParameters[2] = heliTrack->get_R();
 
356
        evEvent.scifiDSTrackerHelicalTrackParameters[3] = heliTrack->get_dsdz();
 
357
        evEvent.scifiDSTrackerHelicalTrackParameters[4] = heliTrack->get_line_sz_c();
 
358
        evEvent.scifiDSTrackerHelicalTrackParameters[5] = startPoint.x();
 
359
        evEvent.scifiDSTrackerHelicalTrackParameters[6] = startPoint.y();
 
360
        evEvent.scifiDSTrackerHelicalTrackParameters[7] = startPoint.z();
 
361
      }
 
362
  }
 
363
}
 
364
 
 
365
int EventViewer::EVExporter::IsVectorNaN(MAUS::ThreeVector v) {
 
366
    if (isnan(v.x()) || isnan(v.y()) || isnan(v.z()))
 
367
        return 1;
 
368
    else
 
369
        return 0;
 
370
}
 
371
 
 
372
} // ~namespace EventViewer