RdScintSignalReconstructor.cc
Go to the documentation of this file.
1 #include "RdScintSignalReconstructor.h"
2 
3 #include <fwk/CentralConfig.h>
4 
5 #include <utl/ErrorLogger.h>
6 #include <utl/Reader.h>
7 #include <utl/config.h>
8 #include <utl/Trace.h>
9 #include <utl/TraceAlgorithm.h>
10 #include <utl/TimeStamp.h>
11 #include <utl/TimeInterval.h>
12 #include <utl/AugerUnits.h>
13 #include <utl/AugerException.h>
14 #include <utl/Math.h>
15 
16 #include <evt/Event.h>
17 #include <evt/ShowerRecData.h>
18 #include <evt/ShowerRRecData.h>
19 #include <evt/ShowerRRecDataQuantities.h>
20 #include <revt/REvent.h>
21 #include <revt/Header.h>
22 #include <revt/Station.h>
23 #include <revt/StationRecData.h>
24 #include <revt/StationRRecDataQuantities.h>
25 #include <revt/StationHeader.h>
26 #include <revt/StationTriggerData.h>
27 #include <revt/Channel.h>
28 #include <revt/ChannelRecData.h>
29 #include <revt/ChannelRRecDataQuantities.h>
30 
31 #include <det/Detector.h>
32 #include <rdet/RDetector.h>
33 #include <rdet/Station.h>
34 #include <rdet/Channel.h>
35 
36 #include <TVector3.h>
37 
38 #define OUT(x) if ((x) <= fInfoLevel) cerr << " "
39 
40 using namespace revt;
41 using namespace evt;
42 using namespace fwk;
43 using namespace utl;
44 using std::cerr;
45 using std::fstream;
46 using std::endl;
47 using std::min;
48 using std::max;
49 
50 
51 namespace RdScintSignalReconstructor {
52 
53  RdScintSignalReconstructor::RdScintSignalReconstructor() :
54  fScintSignalSearchWindowStart(0),
55  fScintSignalSearchWindowStop(0),
56  fNoiseWindowStart(0),
57  fSamplesToCalculateOffset(0),
58  fMinNumberOfScintTop(0),
59  fMinNumberOfScintBottom(0),
60  fMinNumberOfScint(0),
61  fMinSignal(0),
62  fWeightedBaryValues(true),
63  fSimSingleMuonEnergyDepositTop(0),
64  fSimSingleMuonEnergyDepositTopError(0),
65  fSimSingleMuonEnergyDepositBottom(0),
66  fSimSingleMuonEnergyDepositBottomError(0),
67  fInfoLevel(0)
68  { }
69 
70 
71  RdScintSignalReconstructor::~RdScintSignalReconstructor()
72  { }
73 
74 
75  VModule::ResultFlag
76  RdScintSignalReconstructor::Init()
77  {
78  INFO("RdScintSignalReconstructor::Init()");
79 
80  // Read in the configurations of the xml file
81  Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdScintSignalReconstructor");
82  topBranch.GetChild("ScintSignalSearchWindowStart").GetData(fScintSignalSearchWindowStart);
83  topBranch.GetChild("ScintSignalSearchWindowStop").GetData(fScintSignalSearchWindowStop);
84  topBranch.GetChild("NoiseWindowStart").GetData(fNoiseWindowStart);
85  topBranch.GetChild("SamplesToCalculateOffset").GetData(fSamplesToCalculateOffset);
86  topBranch.GetChild("MinSignal").GetData(fMinSignal);
87  topBranch.GetChild("ExcludeStations").GetData(fExcludedStationsName);
88  topBranch.GetChild("MinNumberOfScint").GetData(fMinNumberOfScint);
89  topBranch.GetChild("WeightedBaryValues").GetData(fWeightedBaryValues);
90  topBranch.GetChild("SimSingleMuonEnergyDepositTop").GetData(fSimSingleMuonEnergyDepositTop);
91  topBranch.GetChild("SimSingleMuonEnergyDepositBottom").GetData(fSimSingleMuonEnergyDepositBottom);
92  topBranch.GetChild("SimSingleMuonEnergyDepositTopError").GetData(fSimSingleMuonEnergyDepositTopError);
93  topBranch.GetChild("SimSingleMuonEnergyDepositBottomError").GetData(fSimSingleMuonEnergyDepositBottomError);
94  topBranch.GetChild("FilenameMuonCalibrationTop").GetData(fFilenameMuonCalibrationTop);
95  topBranch.GetChild("FilenameMuonCalibrationBottom").GetData(fFilenameMuonCalibrationBottom);
96 
97  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
98 
99  // read in scintillator calibrations
100  OUT(eDebug) << "read in scintillator calibration top" << endl;
101  fstream fin;
102  char line[256+1];
103  std::string filename = fFilenameMuonCalibrationTop;
104  fin.open(filename.c_str(), std::ios::in); // opening data
105  fin.getline(line, 256); // read/skip first line
106 
107  OUT(eDebug) << "loop through file... filename = " << filename << endl;
108  if (!fin.good()) {
109  WARNING("File TOP_Mp.txt could not be opened");
110  return eFailure;
111  }
112  while (fin.good()) {
113  std::vector<double> vTmp;
114  int stationId;
115  double signalMuonPeak, signalMuonPeakError;
116  fin >> stationId >> signalMuonPeak >> signalMuonPeakError;
117  fin.getline(line, 256); // read rest of the lines
118  stationId += 100; // convert to Offline ids
119  CalibrationData calib;
120  calib.signalMuonPeak = signalMuonPeak;
121  calib.signalMuonPeakError = signalMuonPeakError;
122  fCalibrationDataTop.insert(std::pair<int, CalibrationData> (stationId, calib));
123  OUT(eDebug) << stationId << "\t" << signalMuonPeak << "\t" << signalMuonPeakError << endl;
124  }
125  fin.close();
126 
127  // read in scintillator calibrations
128  OUT(eDebug) << "read in scintillator calibration bottom" << endl;
129  filename = fFilenameMuonCalibrationBottom;
130  fin.open(filename.c_str(), std::ios::in); // opening data
131  fin.getline(line, 256); // read/skip first line
132 
133  OUT(eDebug) << "loop through file... filename = " << filename << endl;
134  if (!fin.good()) {
135  WARNING("File BOTTOM_Mp.txt could not be opened");
136  return eFailure;
137  }
138  while (fin.good()) {
139  std::vector<double> vTmp;
140  int stationId;
141  double signalMuonPeak, signalMuonPeakError;
142  fin >> stationId >> signalMuonPeak >> signalMuonPeakError;
143  fin.getline(line, 256); // read rest of the lines
144  stationId += 100; // convert to Offline ids
145  CalibrationData calib;
146  calib.signalMuonPeak = signalMuonPeak;
147  calib.signalMuonPeakError = signalMuonPeakError;
148  fCalibrationDataBottom.insert(std::pair<int, CalibrationData> (stationId, calib));
149  OUT(eDebug) << stationId << "\t" << signalMuonPeak << "\t" << signalMuonPeakError << endl;
150  }
151  fin.close();
152 
153  return eSuccess;
154  }
155 
156 
157  VModule::ResultFlag
158  RdScintSignalReconstructor::Run(evt::Event& event)
159  {
160  //if (fInfoLevel >= eIntermediate) {
161  INFO("RdScintSignalReconstructor::Run()");
162  //}
163 
164  // Check if there are events at all
165  if (!event.HasREvent()) {
166  WARNING("No radio event found!");
167  return eContinueLoop;
168  }
169 
170  std::stringstream fMessage;
171 
172  REvent& rEvent = event.GetREvent();
173  //TimeStamp eventTime = rEvent.GetHeader().GetTime();
174 
175  unsigned int numberOfTopScintWithPulseFound = 0;
176  unsigned int numberOfBottomScintWithPulseFound = 0;
177  unsigned int numberOfScintStationsWithPulseFound = 0;
178 
179  // loop through stations and for every station through every channel
180  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
181  Station& station = *sIt;
182  int stationId = station.GetId();
183 
184  // RRecStation event should have been generated by RdEventInitializer
185  if (!station.HasRecData()) {
186  WARNING("Station has no RecData! Please call RdEventInitializer first!");
187  return eFailure;
188  }
189  StationRecData& recStation = station.GetRecData();
190 
191  // time that is added to all times determined in a trace in ns
192  double relTime = recStation.GetParameter(eTraceStartTime);
193 
194  std::ostringstream info;
195  info << "Apply pulse finder on station " << station.GetId();
196  OUT(eObscure) << info.str() << endl;
197 
198  // First step, we check if the station belongs to the list of stations to exclude
199  bool ignorestation = false;
200  for (std::vector<std::string>::iterator iditer = fExcludedStationsName.begin();
201  iditer != fExcludedStationsName.end(); ++iditer) {
202  if (rEvent.HasStation(*iditer) &&
203  rEvent.GetStationByName(*iditer).GetId() == sIt->GetId()) {
204  ignorestation = true;
205  break;
206  }
207  }
208  if (ignorestation)
209  continue;
210 
211  double PeakTop = 0;
212  double PeakBottom = 0;
213 
214  //double PeakTopError = 0;
215  //double PeakBottomError = 0;
216 
217  double PeakTopTime = 0;
218  double PeakBottomTime = 0;
219 
220  for (Station::ChannelIterator cIt = station.ChannelsBegin();cIt != station.ChannelsEnd(); ++cIt) {
221  Channel& channel = *cIt;
222 
223  const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
224  const rdet::Channel& detChannel =
225  rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());
226 
227  if (detChannel.GetChannelType() != "ScintillatorTop" &&
228  detChannel.GetChannelType() != "ScintillatorBottom")
229  continue;
230 
231  revt::ChannelTimeSeries& channeltrace = channel.GetChannelTimeSeries();
232  double TraceStop = (channeltrace.GetSize() - 1) * channeltrace.GetBinning();
233 
234  double NoiseWindowStart = fNoiseWindowStart;
235  double ScintSearchWindowStart = fScintSignalSearchWindowStart;
236  double ScintSearchWindowStop = fScintSignalSearchWindowStop;
237  unsigned int samples_offset = fSamplesToCalculateOffset;
238 
239  unsigned int sample = 0;
240  unsigned int NoiseWindowSize = 0;
241 
242  double PeakAmplitude = 0;
243  double IntegratedSignal = 0;
244  double SignalError = 0;
245  double PeakTime = 0;
246  double PeakTimeError = 0;
247  double RMSNoise = 0;
248 
249  // Work on the trace of this channel
250  if (!station.HasChannel(channel.GetId()))
251  continue;
252 
253  // clip windows at trace boundaries
254  ScintSearchWindowStart = min(max(0.0, ScintSearchWindowStart), TraceStop);
255  ScintSearchWindowStop = max(min(ScintSearchWindowStop, TraceStop),ScintSearchWindowStart);
256  //windowfinder
257  double SignalWindowStart = ScintSearchWindowStart;
258  double SignalWindowStop = ScintSearchWindowStop;
259 
260  //pulsefinder
261  Pulsefinder(channeltrace, PeakAmplitude, PeakTime, PeakTimeError, SignalWindowStart, SignalWindowStop, sample, samples_offset);
262 
263  //check if signal is big enough
264  if (PeakAmplitude < fMinSignal || sample == 0)
265  continue;
266 
267  //signalwindowfinder (between zero crossings in front of and after the peak
268  Signalwindowfinder(channeltrace, ScintSearchWindowStart, ScintSearchWindowStop, SignalWindowStart, SignalWindowStop, sample, samples_offset);
269 
270  double SignalIntegrationWindowSize = (SignalWindowStop - SignalWindowStart);
271 
272  // calculate noise
273  NoiseWindowStart = min(max(0.0, NoiseWindowStart), TraceStop);
274  NoiseWindowSize = SignalIntegrationWindowSize;
275 
276  //calculate noise integral
277  Noisefinder(channeltrace, NoiseWindowStart, NoiseWindowSize, SignalError, samples_offset, RMSNoise);
278 
279  PulseFixedWindowIntegrator(channeltrace, SignalIntegrationWindowSize, IntegratedSignal, SignalWindowStart, SignalWindowStop, samples_offset);
280 
281  IntegratedSignal = IntegratedSignal * channeltrace.GetBinning();
282 
283  SignalError = SignalError * channeltrace.GetBinning();
284 
285  if (!channel.HasRecData())
286  channel.MakeRecData();
287 
288  ChannelRecData& recChannel = channel.GetRecData();
289 
290  if (detChannel.GetChannelType() == "ScintillatorTop") {
291  double integratedSignalToVEMTop = 0;
292  if (fCalibrationDataTop.find(stationId) != fCalibrationDataTop.end()) {
293  // found
294  OUT(eDebug) << "Station No. " << stationId << endl;
295  integratedSignalToVEMTop = fCalibrationDataTop[stationId].signalMuonPeak;
296  }
297  if (integratedSignalToVEMTop == 0) {
298  integratedSignalToVEMTop = fCalibrationDataTop[10099].signalMuonPeak;
299  }
300  double integratedSignalToVEMTopError =
301  fCalibrationDataTop[stationId].signalMuonPeakError;
302  double integratedSignalToDepEnergyTop =
303  fSimSingleMuonEnergyDepositTop / integratedSignalToVEMTop;
304  double integratedSignalToDepEnergyTopError =
305  sqrt(pow(fSimSingleMuonEnergyDepositTopError/fSimSingleMuonEnergyDepositTop, 2) +
306  pow(integratedSignalToVEMTopError/integratedSignalToVEMTop, 2)) *
307  integratedSignalToDepEnergyTop;
308 
309  ++numberOfTopScintWithPulseFound;
310  PeakTop = PeakAmplitude;
311  //PeakTopError = RMSNoise;
312  PeakTopTime = relTime + PeakTime;
313 
314  //for now only Top Scintillators VEM and deposited Energy are used later on, therefore fill Station here
315 
316  recStation.SetParameter(eAERAScintStationVEM,
317  IntegratedSignal / integratedSignalToVEMTop);
318  //give VEM unscertainty as sqrt(eAERAScintVEM) for now... should be modified at some point
319  //recStation.SetParameterError(eAERAScintVEM,sqrt(pow(SignalError/IntegratedSignal, 2) + pow(fIntToVEMerror/fIntToVEMTop, 2)*IntegratedSignal/fIntToVEMTop);
320  recStation.SetParameterError(eAERAScintStationVEM,
321  sqrt(IntegratedSignal / integratedSignalToVEMTop));
322  recStation.SetParameter(eAERAScintStationDepositedEnergy,
323  IntegratedSignal * integratedSignalToDepEnergyTop);
324  recStation.SetParameterError(
325  eAERAScintStationDepositedEnergy,
326  sqrt(
327  pow(SignalError / IntegratedSignal, 2) +
328  pow(integratedSignalToDepEnergyTopError / integratedSignalToDepEnergyTop, 2)) *
329  IntegratedSignal * integratedSignalToDepEnergyTop
330  );
331 
332  recChannel.SetParameter(eAERAScintVEM, IntegratedSignal / integratedSignalToVEMTop);
333  OUT(eDebug) << "Integrated Top Signal " << IntegratedSignal << " Top Signal in VEM " << IntegratedSignal / integratedSignalToVEMTop << " Calibration " << integratedSignalToVEMTop << endl;
334  //give VEM unscertainty as sqrt(eAERAScintVEM) for now... should be modified at some point
335  // recChannel.SetParameterError(eAERAScintVEM,sqrt(pow(SignalError/IntegratedSignal, 2) + pow(fIntToVEMerror/fIntToVEMTop, 2)*IntegratedSignal/fIntToVEMTop);
336  recChannel.SetParameterError(eAERAScintVEM,
337  sqrt(IntegratedSignal / integratedSignalToVEMTop));
338 
339  recChannel.SetParameter(eAERAScintDepositedEnergy,
340  IntegratedSignal * integratedSignalToDepEnergyTop);
341  recChannel.SetParameterError(
342  eAERAScintDepositedEnergy,
343  sqrt(
344  pow(SignalError / IntegratedSignal, 2) +
345  pow(integratedSignalToDepEnergyTopError / integratedSignalToDepEnergyTop, 2)) *
346  IntegratedSignal * integratedSignalToDepEnergyTop
347  );
348  }
349  if (detChannel.GetChannelType() == "ScintillatorBottom") {
350  double integratedSignalToVEMBottom = 0;
351  if (fCalibrationDataTop.find(stationId) != fCalibrationDataTop.end()) {
352  //found
353  integratedSignalToVEMBottom = fCalibrationDataBottom[stationId].signalMuonPeak;
354  }
355  if (integratedSignalToVEMBottom == 0) {
356  integratedSignalToVEMBottom = fCalibrationDataBottom[10099].signalMuonPeak;
357  }
358  double integratedSignalToVEMBottomError =
359  fCalibrationDataBottom[stationId].signalMuonPeakError;
360  double integratedSignalToDepEnergyBottom =
361  fSimSingleMuonEnergyDepositBottom / integratedSignalToVEMBottom;
362  double integratedSignalToDepEnergyBottomError =
363  sqrt(
364  pow(fSimSingleMuonEnergyDepositBottomError/fSimSingleMuonEnergyDepositBottom, 2) +
365  pow(integratedSignalToVEMBottomError/integratedSignalToVEMBottom, 2)) *
366  integratedSignalToDepEnergyBottom;
367  ++numberOfBottomScintWithPulseFound;
368  PeakBottom = PeakAmplitude;
369  //PeakBottomError = RMSNoise;
370  PeakBottomTime = relTime + PeakTime;
371 
372  recChannel.SetParameter(eAERAScintVEM,(IntegratedSignal/integratedSignalToVEMBottom));
373  OUT(eDebug) << "Integrated Bottom Signal " << IntegratedSignal << " Bottom Signal in VEM " << IntegratedSignal / integratedSignalToVEMBottom << endl;
374  //give VEM unscertainty as sqrt(eAERAScintVEM) for now... should be modified at some point
375  //recChannel.SetParameterError(eAERAScintVEM,sqrt(pow(SignalError/IntegratedSignal, 2) + pow(fIntToVEMerror/fIntToVEMBottom, 2)*IntegratedSignal/fIntToVEMBottom);
376  recChannel.SetParameterError(eAERAScintVEM,(sqrt(IntegratedSignal/integratedSignalToVEMBottom)));
377  recChannel.SetParameter(eAERAScintDepositedEnergy,IntegratedSignal*integratedSignalToDepEnergyBottom);
378  recChannel.SetParameterError(eAERAScintDepositedEnergy,sqrt(pow(SignalError/IntegratedSignal, 2) + pow(integratedSignalToDepEnergyBottomError/integratedSignalToDepEnergyBottom, 2))*IntegratedSignal*integratedSignalToDepEnergyBottom);
379 
380  }
381 
382  recChannel.SetParameter(eAERAScintSignalTime, relTime + PeakTime);
383  recChannel.SetParameterError(eAERAScintSignalTime, PeakTimeError);
384 
385  recChannel.SetParameter(eAERAScintNoise, RMSNoise);
386 
387  recChannel.SetParameter(eAERAScintIntergratedSignal, IntegratedSignal);
388  recChannel.SetParameterError(eAERAScintIntergratedSignal, SignalError);
389 
390  recChannel.SetParameter(eAERAScintSignalHeight, PeakAmplitude);
391  recChannel.SetParameterError(eAERAScintSignalHeight, RMSNoise);
392 
393  recChannel.SetParameter(eAERAScintSignalWindowStart, SignalWindowStart);
394  recChannel.SetParameter(eAERAScintSignalWindowStop, SignalWindowStop);
395 
396  recChannel.SetParameter(eAERAScintMinSignal, fMinSignal);
397 
398  }
399 
400  if (PeakTop != 0 && PeakBottom != 0) {
401 
402  if (max(PeakTop,PeakBottom) == PeakTop) {
403 
404  recStation.SetParameter(eAERAScintStationSignalTime, PeakTopTime);
405  recStation.SetParameter(eAERAScintStationSignalHeight, PeakTop);
406  recStation.SetParameterError(eAERAScintStationSignalTime, 5);
407 
408  } else {
409  recStation.SetParameter(eAERAScintStationSignalTime, PeakBottomTime);
410  recStation.SetParameter(eAERAScintStationSignalHeight, PeakBottom);
411  recStation.SetParameterError(eAERAScintStationSignalTime, 5);
412  }
413 
414  ++numberOfScintStationsWithPulseFound;
415 
416  } else if (PeakTop != 0) {
417 
418  recStation.SetParameter(eAERAScintStationSignalTime, PeakTopTime);
419  recStation.SetParameter(eAERAScintStationSignalHeight, PeakTop);
420  recStation.SetParameterError(eAERAScintStationSignalTime, 5);
421 
422  ++numberOfScintStationsWithPulseFound;
423 
424  } else if (PeakBottom != 0) {
425 
426  recStation.SetParameter(eAERAScintStationSignalTime, PeakBottomTime);
427  recStation.SetParameter(eAERAScintStationSignalHeight, PeakBottom);
428  recStation.SetParameterError(eAERAScintStationSignalTime, 5);
429 
430  ++numberOfScintStationsWithPulseFound;
431  }
432  }
433 
434  if (!event.HasRecShower())
435  event.MakeRecShower();
436  evt::ShowerRecData& shower = event.GetRecShower();
437 
438  if (!shower.HasRRecShower())
439  shower.MakeRRecShower();
440  evt::ShowerRRecData& rShower = shower.GetRRecShower();
441 
442  rShower.SetParameter(eNumberOfAERATopScintWithPulseFound, numberOfTopScintWithPulseFound);
443  rShower.SetParameter(eNumberOfAERABottomScintWithPulseFound, numberOfBottomScintWithPulseFound);
444  rShower.SetParameter(eNumberOfAERAScintStationWithPulseFound, numberOfScintStationsWithPulseFound);
445 
446  //compute the barycenter of the Scint similar to RdEventInitializer and RdPlaneWaveFit
447  //compute the barytime of the Scint similar to RdEventInitializer and RdPlaneWaveFit
448 
449  const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
450 
451  if (numberOfScintStationsWithPulseFound >= fMinNumberOfScint) {
452 
453  utl::Point BaryCenter;
454  ComputeBaryCenter(rEvent, BaryCenter);
455 
456  rShower.SetParameter(eBaryCenterAERAScintX, BaryCenter.GetX(referenceCS));
457  rShower.SetParameter(eBaryCenterAERAScintY, BaryCenter.GetY(referenceCS));
458  rShower.SetParameter(eBaryCenterAERAScintZ, BaryCenter.GetZ(referenceCS));
459 
460  double BaryTime = 0;
461  ComputeBaryTime(rEvent, BaryTime);
462 
463  rShower.SetParameter(eBaryTimeAERAScint, BaryTime);
464  }
465 
466  OUT(eFinal)
467  << " Number of Stations with ScintSignal = " << numberOfScintStationsWithPulseFound
468  << endl;
469 
470  return eSuccess;
471 
472  }
473 
474 
475  VModule::ResultFlag
476  RdScintSignalReconstructor::Finish()
477  {
478  INFO("RdScintSignalReconstructor::Finish()");
479 
480  return eSuccess;
481  }
482 
483 
484  void
485  RdScintSignalReconstructor::Signalwindowfinder(const revt::ChannelTimeSeries channeltrace,
486  double ScintSearchWindowStart,
487  double ScintSearchWindowStop,
488  double& SignalWindowStart,
489  double& SignalWindowStop,
490  unsigned int SignalTime,
491  unsigned int samples_offset)
492  const
493  {
494  unsigned int startSignalwindowfinder = ScintSearchWindowStart/ channeltrace.GetBinning();
495  unsigned int stopSignalwindowfinder = ScintSearchWindowStop/ channeltrace.GetBinning();
496 
497  double offset = TraceAlgorithm::Mean(channeltrace,startSignalwindowfinder-samples_offset, startSignalwindowfinder);
498 
499  SignalWindowStart = SignalTime;
500  SignalWindowStop = SignalTime;
501 
502  while (SignalWindowStart > startSignalwindowfinder) {
503  if ((channeltrace[SignalWindowStart]-offset)*(-1.) >= 0.) {
504  --SignalWindowStart;
505  } else if ((channeltrace[SignalWindowStart]-offset)*(-1.) < 0. && SignalWindowStart == SignalTime) {
506  break;
507  } else {
508  ++SignalWindowStart;
509  break;
510  }
511  }
512 
513  while (SignalWindowStop <= stopSignalwindowfinder) {
514  if ((channeltrace[SignalWindowStop]-offset)*(-1.) >= 0.) {
515  ++SignalWindowStop;
516  } else if ((channeltrace[SignalWindowStop]-offset)*(-1.) < 0. && SignalWindowStop == SignalTime) {
517  break;
518  } else {
519  --SignalWindowStop;
520  break;
521  }
522  }
523 
524  SignalWindowStart = SignalWindowStart*channeltrace.GetBinning();
525  SignalWindowStop = SignalWindowStop*channeltrace.GetBinning();
526  }
527 
528 
529  void
530  RdScintSignalReconstructor::Noisefinder(const revt::ChannelTimeSeries& channeltrace,
531  double NoiseWindowStart, double NoiseWindowSize,
532  double& SignalError, unsigned int samples_offset,
533  double& RMSNoise)
534  const
535  {
536  unsigned int start = NoiseWindowStart/channeltrace.GetBinning();
537  unsigned int size = NoiseWindowSize/channeltrace.GetBinning();
538  RMSNoise = TraceAlgorithm::RootMeanSquare(channeltrace, start, start+size);
539 
540  double offset = TraceAlgorithm::Mean(channeltrace, start-samples_offset, start);
541 
542  for (unsigned int i = start; i <= start+size; ++i) {
543  SignalError = SignalError + ((channeltrace[i]-offset)*(-1.));
544  }
545  }
546 
547 
548  void
549  RdScintSignalReconstructor::Pulsefinder(const revt::ChannelTimeSeries& channeltrace,
550  double& PeakAmplitude, double& PeakTime, double& PeakTimeError,
551  double SignalWindowStart, double SignalWindowStop,
552  unsigned int& sample, unsigned int samples_offset)
553  const
554  {
555  PeakAmplitude = 0.0;
556  unsigned int start = SignalWindowStart/channeltrace.GetBinning();
557  unsigned int stop = SignalWindowStop/channeltrace.GetBinning();
558 
559  double offset = TraceAlgorithm::Mean(channeltrace,start-samples_offset, start);
560  for (unsigned int i = start; i <= stop; i++) {
561  if ((channeltrace[i]-offset)*(-1.) > PeakAmplitude) {
562  PeakAmplitude = (channeltrace[i]-offset)*(-1.);
563  PeakTime = i * channeltrace.GetBinning();
564  sample = i;
565  }
566  }
567 
568 //#warning Error of PeakTime is currently set to bin size
569  PeakTimeError = channeltrace.GetBinning(); // Binsize as error, has to be changed later
570  }
571 
572 
573  void
574  RdScintSignalReconstructor::PulseFixedWindowIntegrator(const revt::ChannelTimeSeries& channeltrace,
575  double IntegrationTime, double& IntegratedSignal,
576  double SignalWindowStart, double SignalWindowStop,
577  unsigned int samples_offset)
578  const
579  {
580  IntegratedSignal = 0.;
581 
582  unsigned int start = SignalWindowStart/channeltrace.GetBinning();
583  unsigned int stop = SignalWindowStop/channeltrace.GetBinning();
584 
585  double offset = TraceAlgorithm::Mean(channeltrace, start-samples_offset, start);
586 
587  for (unsigned int i = start; i <= stop; i++) {
588  IntegratedSignal += (channeltrace[i]-offset)*(-1);
589  }
590 
591  OUT(eObscure) << "RdScintSignalReconstructor::PulseFixedWindowIntegrator():" << endl;
592  OUT(eObscure) << " IntegrationTime = " << IntegrationTime << " ns" << endl;
593  OUT(eObscure) << " SignalWindowStart = " << SignalWindowStart/ns << " ns" << endl;
594  OUT(eObscure) << " SignalWindowStop = " << SignalWindowStop/ns << " ns" << endl;
595  OUT(eObscure) << " IntegratedSignal = " << IntegratedSignal/(micro*volt/meter)/(micro*volt/meter) << "(\\mu V)^2" << endl;
596  }
597 
598 
599  void
600  RdScintSignalReconstructor::ComputeBaryCenter(const revt::REvent& rEvent,
601  utl::Point& BaryCenter)
602  const
603  {
604  const det::Detector& detector = det::Detector::GetInstance();
605  const rdet::RDetector& rDetector = detector.GetRDetector();
606  const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();
607 
608  // the absolute points in time and space
609  const Point siteOrigin(0, 0, 0, siteCS);
610 
611  Vector barySum(0, 0, 0, siteCS);
612 
613  double weightSum = 0;
614  unsigned int nstations = 0;
615 
616  for (REvent::ConstStationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
617 
618  const StationRecData& sRec = sIt->GetRecData();
619  const rdet::Station& dStation = rDetector.GetStation(*sIt);
620 
621  double signal = 0;
622 
623  if (!sRec.HasParameter(eAERAScintStationSignalHeight)) {
624  continue;
625  }
626 
627  signal = sRec.GetParameter(eAERAScintStationSignalHeight);
628 
629  double weight = sqrt(signal);
630  weightSum += weight;
631  if (fWeightedBaryValues)
632  barySum += weight * (dStation.GetPosition() - siteOrigin);
633  else
634  barySum += (dStation.GetPosition() - siteOrigin);
635  ++nstations;
636 
637  }
638 
639  if (fWeightedBaryValues)
640  barySum /= weightSum;
641  else
642  barySum /= nstations;
643  BaryCenter = barySum + siteOrigin;
644  }
645 
646 
647  void
648  RdScintSignalReconstructor::ComputeBaryTime(const revt::REvent& rEvent, double& BaryTime)
649  const
650  {
651  double weightSum = 0;
652  double TimeSum = 0;
653  double signal = 0;
654  unsigned int nstations = 0;
655  int time = 0;
656 
657  for (REvent::ConstStationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
658 
659  const StationRecData& sRec = sIt->GetRecData();
660 
661  if (!sRec.HasParameter(eAERAScintStationSignalHeight)) {
662  continue;
663  }
664 
665  signal = sRec.GetParameter(eAERAScintStationSignalHeight);
666 
667  time = sRec.GetParameter(eAERAScintStationSignalTime);
668 
669  double weight = sqrt(signal);
670  weightSum += weight;
671  if (fWeightedBaryValues)
672  TimeSum += weight * time;
673  else
674  TimeSum += time;
675  ++nstations;
676 
677  }
678 
679  if (fWeightedBaryValues)
680  TimeSum /= weightSum;
681  else
682  TimeSum /= nstations;
683  BaryTime = TimeSum;
684 
685  }
686 
687 }
revt::ChannelRecData
Class to access channel level reconstructed data.
Definition: REvent/ChannelRecData.h:23
utl::Branch::GetTopBranch
Branch GetTopBranch() const
Definition: Branch.cc:63
revt::Channel::MakeRecData
void MakeRecData()
Make channel reconstructed data object.
Definition: REvent/Channel.cc:10
revt::StationRecData
Class to access station level reconstructed data.
Definition: REvent/StationRecData.h:26
evt::ShowerRRecData::SetParameter
void SetParameter(Parameter i, double value, bool lock=true)
Definition: ShowerRRecData.h:83
revt::StationRecData::SetParameterError
void SetParameterError(Parameter i, double value, bool lock=true)
Definition: REvent/StationRecData.h:62
utl::Point
Point object.
Definition: Point.h:32
revt::Channel::GetId
int GetId() const
Return Id of the Channel.
Definition: REvent/Channel.h:31
fwk::VModule::eSuccess
Report success to RunController.
Definition: VModule.h:62
rdet::Station
Detector description interface for Station-related data.
Definition: RDetector/Station.h:50
RdScintSignalReconstructor::RdScintSignalReconstructor::Noisefinder
void Noisefinder(const revt::ChannelTimeSeries &channeltrace, double NoiseWindowStart, double NoiseWindowSize, double &SignalError, unsigned int samples_offset, double &RMSNoise) const
Definition: RdScintSignalReconstructor.cc:530
revt::Station::GetRecData
StationRecData & GetRecData()
Get station level reconstructed data.
Definition: REvent/Station.h:95
evt::ShowerRecData
Interface class to access Shower Reconstructed parameters.
Definition: ShowerRecData.h:33
revt::REvent::GetStationByName
Station & GetStationByName(const std::string &name)
retrieve station by name, throw utl::NonExistentComponentException if n.a.
Definition: REvent.h:194
evt::Event::HasRecShower
bool HasRecShower() const
Definition: Framework/Event/Event.h:55
revt::REvent
Interface class to access to the Radio part of an event.
Definition: REvent.h:42
fwk::VModule::eContinueLoop
Skip remaining modules in the current loop and continue with next iteration of the loop...
Definition: VModule.h:68
evt::ShowerRRecData
Interface class to access to the RD Reconstruction of a Shower.
Definition: ShowerRRecData.h:32
revt::ChannelRecData::SetParameter
void SetParameter(Parameter i, double value, bool lock=true)
Definition: REvent/ChannelRecData.h:53
utl::Trace::GetBinning
double GetBinning() const
size of one slot
Definition: Trace.h:138
galactic::meter
const double meter
Definition: GalacticUnits.h:29
revt::Channel::GetRecData
ChannelRecData & GetRecData()
Get channel level reconstructed data.
Definition: REvent/Channel.h:90
INFO
#define INFO(message)
Macro for logging informational messages.
Definition: ErrorLogger.h:161
RdScintSignalReconstructor::RdScintSignalReconstructor::Signalwindowfinder
void Signalwindowfinder(const revt::ChannelTimeSeries channeltrace, double ScintSearchWindowStart, double ScintSearchWindowStop, double &SignalWindowStart, double &SignalWindowStop, unsigned int sample, unsigned int samples_offset) const
Definition: RdScintSignalReconstructor.cc:485
revt::REvent::StationsEnd
StationIterator StationsEnd()
Definition: REvent.h:130
revt::REvent::StationsBegin
StationIterator StationsBegin()
Definition: REvent.h:128
fwk::CoordinateSystemRegistry::Init
void Init()
Initialise the registry.
Definition: CoordinateSystemRegistry.cc:58
rdet::Channel
Detector description interface for Channel-related data.
Definition: RDetector/Channel.h:27
utl::Branch::GetChild
Branch GetChild(const std::string &childName) const
Get child of this Branch by child name.
Definition: Branch.cc:211
det::Detector::GetSiteCoordinateSystem
utl::CoordinateSystemPtr GetSiteCoordinateSystem() const
Get the coordinate system for the site.
Definition: Detector.h:137
revt::Station::ChannelsBegin
ChannelIterator ChannelsBegin()
begin Channel iterator for read/write
Definition: REvent/Station.h:66
std::pow
double pow(const double x, const unsigned int i)
Definition: MuonArrivalTime.h:17
revt::REvent::StationIterator
boost::filter_iterator< StationFilter, AllStationIterator > StationIterator
Iterator over all (non-exculded) stations.
Definition: REvent.h:125
evt::Event::HasREvent
bool HasREvent() const
Definition: Framework/Event/Event.h:70
revt::Channel::HasRecData
bool HasRecData() const
Check whether channel reconstructed data exists.
Definition: REvent/Channel.h:96
evt::ShowerRecData::GetRRecShower
ShowerRRecData & GetRRecShower()
Definition: ShowerRecData.cc:134
revt::Channel::GetChannelTimeSeries
ChannelTimeSeries & GetChannelTimeSeries()
retrieve Channel Time Series (write access, only use this if you intend to change the data) ...
Definition: REvent/Channel.h:53
revt::Station::ChannelsEnd
ChannelIterator ChannelsEnd()
end Channel iterator for read/write
Definition: REvent/Station.h:70
rdet::RDetector
Detector description interface for RDetector-related data.
Definition: RDetector.h:46
revt::Station::HasChannel
bool HasChannel(const int pmtId) const
Check if a particular Channel object exists.
Definition: REvent/Station.cc:100
utl::CoordinateSystemPtr
boost::shared_ptr< const CoordinateTransformer > CoordinateSystemPtr
Shared pointer for coordinate systems.
Definition: CoordinateSystemPtr.h:9
rdet::Channel::GetChannelType
const std::string & GetChannelType() const
Get description of Channel Type.
Definition: RDetector/Channel.cc:82
max
#define max(a, b)
Definition: DLECorrectionGG/DLECorrection.cc:30
utl::Branch
Class representing a document branch.
Definition: Branch.h:107
revt::Station
class to hold data at the radio Station level.
Definition: REvent/Station.h:39
utl::BasicVector::GetX
double GetX(const CoordinateSystemPtr &coordinateSystem) const
Definition: BasicVector.h:206
ns
const double ns
Definition: PhysicsConstants.h:2
revt::ChannelRecData::SetParameterError
void SetParameterError(Parameter i, double value, bool lock=true)
Definition: REvent/ChannelRecData.h:56
evt::ShowerRecData::HasRRecShower
bool HasRRecShower() const
Definition: ShowerRecData.cc:149
det::Detector
Top of the hierarchy of the detector description interface.
Definition: Detector.h:81
utl::Trace::GetSize
SizeType GetSize() const
Definition: Trace.h:156
WARNING
#define WARNING(message)
Macro for logging warning messages.
Definition: ErrorLogger.h:163
utl::Branch::GetData
void GetData(bool &b) const
Overloads of the GetData member template function.
Definition: Branch.cc:644
revt::Station::GetId
int GetId() const
Get the station Id.
Definition: REvent/Station.h:56
OUT
#define OUT(x)
Definition: RdScintSignalReconstructor.cc:38
revt::REvent::ConstStationIterator
boost::filter_iterator< StationFilter, ConstAllStationIterator > ConstStationIterator
Definition: REvent.h:126
revt::StationRecData::GetParameter
double GetParameter(const Parameter i) const
Definition: REvent/StationRecData.h:32
utl::BasicVector::GetY
double GetY(const CoordinateSystemPtr &coordinateSystem) const
Definition: BasicVector.h:209
RdScintSignalReconstructor::RdScintSignalReconstructor::Pulsefinder
void Pulsefinder(const revt::ChannelTimeSeries &channeltrace, double &PeakAmplitude, double &PeakTime, double &PeakTimeError, double SignalWindowStart, double SignalWindowStop, unsigned int &sample, unsigned int samples_offset) const
Definition: RdScintSignalReconstructor.cc:549
fwk::VModule::ResultFlag
ResultFlag
Flag returned by module methods to the RunController.
Definition: VModule.h:60
revt::StationRecData::HasParameter
bool HasParameter(const Parameter i) const
Definition: REvent/StationRecData.h:50
revt::StationRecData::SetParameter
void SetParameter(Parameter i, double value, bool lock=true)
Definition: REvent/StationRecData.h:59
utl::Vector
Vector object.
Definition: Vector.h:30
revt::Channel
Class that holds the data associated to an individual radio channel.
Definition: REvent/Channel.h:27
fwk::VModule::eFailure
Report failure to RunController, causing RunController to terminate execution.
Definition: VModule.h:64
utl::micro
constexpr double micro
Definition: AugerUnits.h:65
det::Detector::GetRDetector
const rdet::RDetector & GetRDetector() const
Definition: Detector.cc:143
RdScintSignalReconstructor::RdScintSignalReconstructor::PulseFixedWindowIntegrator
void PulseFixedWindowIntegrator(const revt::ChannelTimeSeries &channeltrace, double IntegrationTime, double &IntegratedSignal, double SignalWindowStart, double SignalWindowStop, unsigned int samples_offset) const
Definition: RdScintSignalReconstructor.cc:574
revt::REvent::HasStation
bool HasStation(const int stationId) const
Check whether station exists.
Definition: REvent.cc:132
un2::Mean
double Mean(const std::vector< double > &v)
Definition: Functions.h:31
filename
char * filename
Definition: dump1090.h:266
revt::Station::HasRecData
bool HasRecData() const
Check whether station reconstructed data exists.
Definition: REvent/Station.h:100
utl::BasicVector::GetZ
double GetZ(const CoordinateSystemPtr &coordinateSystem) const
Definition: BasicVector.h:212
rdet::Station::GetPosition
utl::Point GetPosition() const
Tank position in Site Cartesian Coordinates.
Definition: RDetector/Station.cc:127
rdet::RDetector::GetStation
const Station & GetStation(const int stationId) const
Get station by Station Id.
Definition: RDetector.cc:141
galactic::volt
const double volt
Definition: GalacticUnits.h:38
RdScintSignalReconstructor::RdScintSignalReconstructor::ComputeBaryTime
void ComputeBaryTime(const revt::REvent &rEvent, double &BaryTime) const
Definition: RdScintSignalReconstructor.cc:648
RdScintSignalReconstructor::RdScintSignalReconstructor::ComputeBaryCenter
void ComputeBaryCenter(const revt::REvent &rEvent, utl::Point &BaryCenter) const
Definition: RdScintSignalReconstructor.cc:600
, generated on Tue Sep 26 2023.