SD2ADST.cc
Go to the documentation of this file.
1 #include <utl/config.h>
2 
3 #include "SD2ADST.h"
4 #include "Config.h"
5 #include "ConversionUtil.h"
6 
7 #include <adst/RecEvent.h>
8 
9 #include <fwk/CentralConfig.h>
10 #include <fwk/LocalCoordinateSystem.h>
11 #include <fwk/CoordinateSystemRegistry.h>
12 
13 #include <det/Detector.h>
14 
15 #include <fdet/FDetector.h>
16 #include <fdet/Pixel.h>
17 #include <fdet/Camera.h>
18 #include <fdet/Eye.h>
19 #include <fdet/Telescope.h>
20 
21 #include <sdet/SDetector.h>
22 #include <sdet/Station.h>
23 #include <sdet/STimeVariance.h>
24 #include <sdet/PMTConstants.h>
25 
26 #include <evt/Event.h>
27 #include <evt/ShowerRecData.h>
28 #include <evt/ShowerSimData.h>
29 #include <evt/ShowerFRecData.h>
30 #include <evt/ShowerSRecData.h>
31 #include <evt/ShowerScintillatorRecData.h>
32 
33 #include <evt/ShowerUnivRecData.h>
34 
35 #include <evt/RiseTime1000.h>
36 #include <sevt/SdFootprintData.h>
37 
38 #include <evt/GaisserHillas4Parameter.h>
39 #include <evt/GaisserHillas6Parameter.h>
40 
41 #include <sevt/SEvent.h>
42 #include <sevt/Header.h>
43 #include <sevt/EventTrigger.h>
44 #include <sevt/Station.h>
45 #include <sevt/StationSimData.h>
46 #include <sevt/PMT.h>
47 #include <sevt/PMTRecData.h>
48 #include <sevt/SortCriteria.h>
49 #include <sevt/Meteo.h>
50 #include <sevt/Scintillator.h>
51 #include <sevt/ScintillatorRecData.h>
52 #include <sevt/ScintillatorSimData.h>
53 
54 #include <utl/Branch.h>
55 #include <utl/ErrorLogger.h>
56 #include <utl/UTCDateTime.h>
57 #include <utl/TimeStamp.h>
58 #include <utl/ModifiedJulianDate.h>
59 #include <utl/MathConstants.h>
60 #include <utl/ReferenceEllipsoid.h>
61 #include <utl/AugerUnits.h>
62 #include <utl/Point.h>
63 #include <utl/Vector.h>
64 #include <utl/UTMPoint.h>
65 #include <utl/Transformation.h>
66 #include <utl/Trace.h>
67 #include <utl/AugerUnits.h>
68 #include <utl/MathConstants.h>
69 #include <utl/PhysicalConstants.h>
70 #include <utl/PhysicalFunctions.h>
71 #include <evt/VGaisserHillasParameter.h>
72 #include <utl/config.h>
73 #include <utl/Particle.h>
74 #include <utl/CoordinateSystemPtr.h>
75 #include <utl/TabulatedFunctionErrors.h>
76 #include <utl/Is.h>
77 
78 #ifdef INCLUDE_HAS
79 # include <sdet/MuonProfile.h>
80 # include <sdet/VTankResponse.h>
81 # include <sdet/TankResponseFactory.h>
82 # include <sdet/EMComponent.h>
83 # include <utl/HASUtilities.h>
84 #endif
85 
86 #include <iostream>
87 #include <string>
88 #include <algorithm>
89 
90 
91 //using namespace evt; // DONT USE NAMESPACE evt !!!!!!
92 using namespace std;
93 using namespace utl;
94 using namespace fwk;
95 using namespace otoa;
96 #ifdef INCLUDE_HAS
97 using namespace HASUtilities;
98 #endif
99 
100 
101 namespace {
102 
103  // square of the perpendicular distance from the axis
104  inline
105  double
106  RPerp2(const Vector& axis, const Vector& station)
107  {
108  const double scal = axis*station;
109  return station*station - scal*scal;
110  }
111 
112 }
113 
114 
115 void
116 SD2ADST::Convert(const evt::Event& event, RecEvent& recEvent)
117  const
118 {
119  FillSEvent(event, recEvent);
120  AddStations(event, recEvent);
121  FillMPD(event, recEvent);
122  FillUniversality(event, recEvent);
123 }
124 
125 
126 void
127 SD2ADST::AddStations(const evt::Event& event, RecEvent& recEvent)
128  const
129 {
130  auto& adstSEvent = recEvent.GetSDEvent();
131  auto& adstSdRecShower = adstSEvent.GetSdRecShower();
132  const auto& cfg = GetConfig();
133 
134  const auto& detector = det::Detector::GetInstance();
135  const auto& sDetector = detector.GetSDetector();
136  const auto& referenceCS = detector.GetReferenceCoordinateSystem();
137  const auto& sevent = event.GetSEvent();
138  sevent.SortStations(sevt::ByDecreasingSignal());
139 
140  const auto& timeVariance = sdet::STimeVariance::GetInstance();
141 
142  for (const auto& station : sevent.StationsRange()) {
143 
144  SdRecStation recStation;
145  recStation.SetId(station.GetId());
146  bool filledSomething = false;
147 
148  const auto& dStation = sDetector.GetStation(station);
149  const auto& sPosition = dStation.GetPosition();
150 
151  if (dStation.IsDense()) {
152  try {
153  recStation.SetIsDense();
154  } catch (const AugerException& e) {
155  ERROR(e.GetMessage());
156  continue;
157  }
158  }
159 
160  // station trigger info
161  if (station.HasTriggerData()) {
162  const auto& trigStation = station.GetTriggerData();
163  if (trigStation.IsTimeOverThreshold())
164  recStation.SetToT();
165  if (trigStation.IsTimeOverThresholdDeconvoluted())
166  recStation.SetToTd();
167  if (trigStation.IsMultiplicityOfPositiveSteps())
168  recStation.SetMoPS();
169  if (trigStation.IsT2Threshold())
170  recStation.SetT2Threshold();
171  if (trigStation.IsT1Threshold())
172  recStation.SetT1Threshold();
173  if (trigStation.IsRDThreshold())
174  recStation.SetRDThreshold();
175  if (trigStation.IsRandom())
176  recStation.SetRandom();
177  recStation.SetPLDVersion(trigStation.GetPLDVersion());
178  recStation.SetPLDTimeOffset(trigStation.GetPLDTimeOffset());
179  recStation.SetT3ErrorCode(trigStation.GetErrorCode());
180  recStation.SetT3Window(trigStation.GetWindowMicroSecond());
181  }
182 
183  if (station.HasCalibData())
184  recStation.SetCalibrationVersion(station.GetCalibData().GetVersion());
185 
186  if (station.IsCandidate())
187  recStation.SetCandidate();
188  else if (station.IsSilent())
189  recStation.SetSilent();
190  else { // only rejected stations remain at this stage
191  recStation.SetRejected();
192  recStation.SetRejectStatus(station.GetRejectionStatus());
193  if (station.HasRecData() && (station.GetRejectionStatus() & eLightning)) {
194  // has data, no silent and rejected --> accidental or lightning or dense
195  adstSEvent.SetLightning();
196  }
197  // no data / silent and rejected --> bad station
198  adstSEvent.AddBadStation(
199  SdBadStation(station.GetId(), station.GetRejectionStatus())
200  );
201  }
202 
203  recStation.SetIsUUB(station.IsUUB());
204 
205  if (station.IsHighGainSaturation())
206  recStation.SetHighGainSaturated();
207  if (station.IsLowGainSaturation())
208  recStation.SetLowGainSaturated();
209 
210  recStation.SetBottomUpResidual(station.GetBottomUpResidual());
211 
212  if (station.HasRecData()) {
213  filledSomething = true;
214 
215  const auto& recData = station.GetRecData();
216 
217  if (recData.GetRecoveredSignal())
218  recStation.SetRecoveredSignal(recData.GetRecoveredSignal(),
219  recData.GetRecoveredSignalError());
220 
221  // set basic station properties
222  recStation.SetTime(recData.GetSignalStartTime().GetGPSSecond(),
223  recData.GetSignalStartTime().GetGPSNanoSecond());
224  const double signal = recData.GetTotalSignal();
225  recStation.SetTotalSignal(signal, recData.GetTotalSignalError());
226  recStation.SetRiseTime(recData.GetRiseTime(), recData.GetRiseTimeRMS());
227  recStation.SetFallTime(recData.GetFallTime(), recData.GetFallTimeRMS());
228  recStation.SetSPDistance(recData.GetSPDistance(), recData.GetSPDistanceError());
229 
230  recStation.SetSignalStartAndEnd(recData.GetSignalStartSlot(), recData.GetSignalEndSlot());
231  recStation.SetShapeParameter(recData.GetShapeParameter(), recData.GetShapeParameterRMS());
232  recStation.SetMuonComponent(recData.GetMuonComponent());
233  recStation.SetTime50(recData.GetT50(), recData.GetT50RMS());
234 
235  // adding scintillator stuff
236  if (station.HasScintillator() && station.GetScintillator().HasRecData()) {
237 
238  recStation.SetHasScintillator(station.HasScintillator());
239 
240  auto& recStSci = recStation.GetScintillator();
241  const auto& sci = station.GetScintillator();
242  const auto& sciRec = sci.GetRecData();
243  recStSci.SetTotalSignal(sciRec.GetTotalSignal(), sciRec.GetTotalSignalError());
244  recStSci.SetRiseTime(sciRec.GetRiseTime(), sciRec.GetRiseTimeError());
245  recStSci.SetLDFResidual(sciRec.GetLDFResidual());
246  recStSci.SetSignalStartAndEnd(sciRec.GetSignalStartSlot(), sciRec.GetSignalEndSlot());
247 
248  if (sci.IsHighGainSaturation())
249  recStSci.SetHighGainSaturated();
250  if (sci.IsLowGainSaturation())
251  recStSci.SetLowGainSaturated();
252 
253  }
254 
255  // finally stuff requiring SdRecShower
256  if (event.HasRecShower() && event.GetRecShower().HasSRecShower()) {
257 
258  const auto& sRecShower = event.GetRecShower().GetSRecShower();
259  const auto& core = sRecShower.GetCorePosition();
260  const auto& coreCS = LocalCoordinateSystem::Create(core);
261  const auto& coreTime = sRecShower.GetCoreTime();
262  const auto& axis = sRecShower.GetAxis();
263  const double curvature = sRecShower.GetCurvature();
264  const double curvRadius = curvature ? 1/curvature : 0;
265  const auto showerCenter = curvRadius * axis;
266  const auto showerOrigin = core + showerCenter;
267 
268  const TimeStamp curvTime0 = coreTime - TimeInterval(curvRadius/kSpeedOfLight);
269 
270  // the time residuals from the plane fit approximation to curvature fit (near core)
271  const auto sRelPos = sPosition - core;
272  const double height = axis * sRelPos;
273  const double radius = (showerCenter - sRelPos).GetMag();
274  const TimeStamp planeFrontTime = coreTime - TimeInterval(height/kSpeedOfLight);
275  const double tCurvature =
276  (curvTime0 + TimeInterval(radius/kSpeedOfLight) - recData.GetSignalStartTime()).GetInterval();
277  const double tPlaneFront = (planeFrontTime - recData.GetSignalStartTime()).GetInterval();
278 
279  if (curvature)
280  recStation.SetCurvatureTimeResidual(tCurvature/nanosecond);
281 
282  recStation.SetPlaneTimeResidual(tPlaneFront/nanosecond);
283  recStation.SetLDFResidual(recData.GetLDFResidual());
284  recStation.SetAzimuthSP(recData.GetAzimuthShowerPlane());
285 
286  const auto stationToOrigin = showerOrigin - sPosition;
287 
288  const double cosThetaStation = stationToOrigin.GetZ(coreCS) / stationToOrigin.GetMag();
289  const double timeVar =
290  timeVariance.GetTimeSigma2(signal, recData.GetT50(), cosThetaStation, sqrt(RPerp2(axis, sRelPos)));
291  recStation.SetTimeVariance(timeVar);
292 
293  }
294 
295  // rise time related things
296  recStation.SetCorrRiseTime(recData.GetCorrectedRiseTime(),
297  recData.GetCorrectedRiseTimeError());
298  recStation.SetRiseTimeRejectCode(recData.GetRiseTimeRejectionCode());
299 
300  bool usedInGlobalMPD = false;
301  // vectors to store the MPDs (one MPD per station)
302  vector<double> mpdSignal;
303  vector<double> mpdDepth;
304  vector<double> mpdSignalEr;
305  vector<double> mpdDepthEr;
306 
307  for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {
308 
309  if (!pmt.HasCalibData())
310  continue;
311  const auto& pmtCalib = pmt.GetCalibData();
312 
313  filledSomething = true;
314  Traces recTraces;
315 
316  recTraces.SetPMTId(pmt.GetId());
317 
318  recTraces.SetDynodeAnodeRatio(pmtCalib.GetGainRatio());
319 
320  {
321  auto gain = sdet::PMTConstants::eHighGain;
322  recTraces.SetBaseline(pmtCalib.GetBaseline(gain), pmtCalib.GetBaselineRMS(gain));
323  gain = sdet::PMTConstants::eLowGain;
324  recTraces.SetBaselineLG(pmtCalib.GetBaseline(gain), pmtCalib.GetBaselineRMS(gain));
325  }
326 
327  if (cfg.StoreCalibHistos()) {
328  const auto& chargeBinning =
329  dStation.GetMuonChargeHistogramBinning<double>(pmt.GetType(), station.GetCalibData().GetVersion());
330  vector<double> meanOfBins;
331  for (unsigned int k = 1, n = chargeBinning.size(); k < n; ++k)
332  meanOfBins.push_back(0.5*(chargeBinning[k] + chargeBinning[k-1]));
333 
334  const auto& charge = pmtCalib.GetMuonChargeHisto();
335 
336  recTraces.SetChargeHistogram(meanOfBins, charge);
337 
338  const auto& peakBinning =
339  dStation.GetMuonPeakHistogramBinning<double>(pmt.GetType(), station.GetCalibData().GetVersion());
340  meanOfBins.clear();
341  for (unsigned int k = 1, n = peakBinning.size(); k < n; ++k)
342  meanOfBins.push_back(0.5*(peakBinning[k] + peakBinning[k-1]));
343 
344  const auto& peak = pmtCalib.GetMuonPeakHisto();
345  recTraces.SetPeakHistogram(meanOfBins, peak);
346  }
347 
348  const double energy = adstSdRecShower.GetEnergy();
349  const bool enoughEnergy = cfg.StoreSDTracesMinEnergy() <= 0 || energy >= cfg.StoreSDTracesMinEnergy();
350  const bool storeVEMTraces = enoughEnergy && cfg.StoreSDTraces() > 0;
351  const bool storeFADCTraces = enoughEnergy && cfg.StoreSDTraces() > 1;
352  const bool storeBaselines = enoughEnergy && cfg.StoreSDTraces() > 2;
353 
354  if (pmt.HasRecData()) {
355  const auto& pmtRec = pmt.GetRecData();
356  if (pmtRec.HasMuonProductionDepth() && energy >= cfg.StoreMPDsMinEnergy()) {
357  const auto& mpds = pmtRec.GetMuonProductionDepth();
358 
359  copy(mpds.YBegin(), mpds.YEnd(), back_inserter(mpdSignal));
360  copy(mpds.XBegin(), mpds.XEnd(), back_inserter(mpdDepth));
361  copy(mpds.YErrBegin(), mpds.YErrEnd(), back_inserter(mpdSignalEr));
362  copy(mpds.XErrBegin(), mpds.XErrEnd(), back_inserter(mpdDepthEr));
363 
364  usedInGlobalMPD = pmtRec.IsUsedInGlobalMPD();
365  }
366 
367  recTraces.SetCharge(pmtRec.GetVEMCharge(), pmtRec.GetVEMChargeError());
368  recTraces.SetOnlineVEMCharge(pmtRec.GetOnlineVEMCharge());
369  recTraces.SetHistogramVEMCharge(pmtRec.GetHistogramVEMCharge());
370  recTraces.SetPeak(pmtRec.GetVEMPeak(), pmtRec.GetVEMPeakError());
371  recTraces.SetOnlineVEMPeak(pmtRec.GetOnlineVEMPeak());
372  recTraces.SetHistogramVEMPeak(pmtRec.GetHistogramVEMPeak());
373  recTraces.SetVEMChargeFromHistogram(pmtRec.IsVEMChargeFromHistogram());
374  recTraces.SetMuonPulseDecayTime(pmtRec.GetMuonPulseDecayTime(),
375  pmtRec.GetMuonPulseDecayTimeError());
376 
377  if (pmtRec.HasVEMTrace()) {
378  recTraces.SetType(GetTraceType(sevt::StationConstants::eTotal));
379  const auto& trace = pmtRec.GetVEMTrace();
380 
381  if (storeVEMTraces) {
382  const vector<float> t(trace.Begin(), trace.End());
383  recTraces.SetVEMComponent(t);
384  }
385 
386  if (pmtRec.GetVEMCharge() > 0) {
387  double xsum = 0;
388  // note: last bin is included in order to be consistent with SdCalibrator
389  for (int i = recData.GetSignalStartSlot(), n = recData.GetSignalEndSlot(); i <= n; ++i)
390  xsum += trace[i];
391  recTraces.SetVEMSignal(xsum * pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge());
392  }
393  }
394  }
395 
396  if (storeFADCTraces && pmt.HasFADCTrace()) {
397 
398  recTraces.SetIsTubeOk(pmt.GetCalibData().IsTubeOk());
399  recTraces.SetIsLowGainOk(pmt.GetCalibData().IsLowGainOk());
400 
401  const auto& traceLG = pmt.GetFADCTrace(sdet::PMTConstants::eLowGain);
402  const vector<short unsigned int> lg(traceLG.Begin(), traceLG.End());
403  recTraces.SetLowGainComponent(lg);
404 
405  const auto& traceHG = pmt.GetFADCTrace(sdet::PMTConstants::eHighGain);
406  vector<short unsigned int> hg(traceHG.Begin(), traceHG.End());
407  recTraces.SetHighGainComponent(hg);
408  }
409 
410  recStation.AddPMTTraces(recTraces);
411  recTraces.Clear();
412 
413  // Store baselines
414  if (storeBaselines && pmt.HasRecData()) {
415  const auto& pmtRec = pmt.GetRecData();
416  if (pmtRec.HasFADCBaseline(sdet::PMTConstants::eLowGain)) {
417  recTraces.SetType(eBaselineLowGain);
418  recTraces.SetPMTId(pmt.GetId());
419  const auto& baselineLG = pmtRec.GetFADCBaseline(sdet::PMTConstants::eLowGain);
420  const vector<float> lg(baselineLG.Begin(), baselineLG.End());
421  recTraces.SetVEMComponent(lg);
422  recStation.AddPMTTraces(recTraces);
423  recTraces.Clear();
424  }
425  if (pmtRec.HasFADCBaseline(sdet::PMTConstants::eHighGain)) {
426  recTraces.SetType(eBaselineHighGain);
427  recTraces.SetPMTId(pmt.GetId());
428  const auto& baselineHG = pmtRec.GetFADCBaseline(sdet::PMTConstants::eHighGain);
429  const vector<float> hg(baselineHG.Begin(), baselineHG.End());
430  recTraces.SetVEMComponent(hg);
431  recStation.AddPMTTraces(recTraces);
432  recTraces.Clear();
433  }
434  }
435 
436  // fill the components traces
437  if (storeVEMTraces && !cfg.StoreMCTraces())
438  continue;
439 
440  if (!pmt.HasRecData())
441  continue;
442  const auto& pmtRec = pmt.GetRecData();
443 
444  for (const auto traces : pmtRec.VEMTracesRange()) {
445  filledSomething = true;
446  const auto sigComp = sevt::StationConstants::SignalComponent(traces.GetLabel());
447  if (sigComp == sevt::StationConstants::eTotal || GetTraceType(sigComp) == eUnknownTrace)
448  continue;
449 
450  recTraces.SetType(GetTraceType(sigComp));
451  recTraces.SetCharge(pmtRec.GetVEMCharge(), pmtRec.GetVEMChargeError());
452  recTraces.SetOnlineVEMCharge(pmtRec.GetOnlineVEMCharge());
453  recTraces.SetHistogramVEMCharge(pmtRec.GetHistogramVEMCharge());
454  recTraces.SetPeak(pmtRec.GetVEMPeak(), pmtRec.GetVEMPeakError());
455  recTraces.SetOnlineVEMPeak(pmtRec.GetOnlineVEMPeak());
456  recTraces.SetHistogramVEMPeak(pmtRec.GetHistogramVEMPeak());
457  recTraces.SetVEMChargeFromHistogram(pmtRec.IsVEMChargeFromHistogram());
458  recTraces.SetPMTId(pmt.GetId());
459  const auto& trace = traces.GetTrace();
460  const vector<float> t(trace.Begin(), trace.End());
461  recTraces.SetVEMComponent(t);
462  recStation.AddPMTTraces(recTraces);
463  recTraces.Clear();
464  }
465  }
466  recStation.SetMPDDepth(mpdDepth, mpdDepthEr);
467  recStation.SetMPDSignal(mpdSignal, mpdSignalEr);
468  recStation.SetUsedInGlobalMPD(usedInGlobalMPD);
469  }
470 
471  if (filledSomething)
472  adstSEvent.AddStation(recStation);
473 
474  // fill the Sd MC station info
475  if (!station.HasSimData())
476  continue;
477 
478  const auto& simStat = station.GetSimData();
479  GenStation myStation;
480 
481  myStation.SetId(station.GetId());
482  const auto& pFront = simStat.GetPlaneFrontTime();
483  myStation.SetPlaneFrontTime(pFront.GetGPSSecond(), pFront.GetGPSNanoSecond());
484  myStation.SetInsideRMinFlag(simStat.IsInsideMinRadius());
485  myStation.SetNumberOfMuons(simStat.GetNumberOfMuons());
486  myStation.SetNumberOfElectrons(simStat.GetNumberOfElectrons());
487  myStation.SetNumberOfPhotons(simStat.GetNumberOfPhotons());
488  myStation.SetTotalParticleCount(simStat.GetTotalParticleCount());
489  myStation.SetMaxNParticles(simStat.GetMaxNParticles());
490  myStation.SetThinning(simStat.GetThinning());
491  myStation.SetThinParticleCount(simStat.GetNParticles());
492  myStation.SetUsedWeight(simStat.GetUsedWeight());
493 
494  if (event.HasSimShower())
495  myStation.SetSPDistance(sPosition.GetRho(event.GetSimShower().GetShowerCoordinateSystem()) / m);
496 
497  if (station.HasScintillator()) {
498  if (station.GetScintillator().HasSimData()) {
499  const auto& sciSim = station.GetScintillator().GetSimData();
500  auto& mySci = myStation.GetScintillator();
501  mySci.SetNumberOfMuons(sciSim.GetNumberOfMuons());
502  mySci.SetNumberOfElectrons(sciSim.GetNumberOfElectrons());
503  mySci.SetNumberOfPhotons(sciSim.GetNumberOfPhotons());
504  mySci.SetTotalParticleCount(sciSim.GetTotalParticleCount());
505  }
506  }
507 
508  if (cfg.StoreSDPETimeDistribution()) {
509 
510  for (const auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {
511 
512  if (pmt.HasSimData()) {
513  if (pmt.GetSimData().HasPETimeDistribution()) {
514  for (const auto& pe : pmt.GetSimData().PETimeDistributionsRange()) {
515 
516  const auto sigComp = sevt::StationConstants::SignalComponent(pe.GetLabel());
517 
518  PETimeDistribution timeDistribution;
519  const auto& timeDist = pe.GetTimeDistribution();
520  for (int timeit = timeDist.GetStart(), n = timeDist.GetStop(); timeit <= n; ++timeit) {
521  if (timeDist[timeit] > 0)
522  timeDistribution.FillMap(timeit, timeDist[timeit]);
523  }
524  timeDistribution.SetPMTId(pmt.GetId());
525  timeDistribution.SetComponent(GetTraceType(sigComp));
526  myStation.AddPETimeDistribution(timeDistribution);
527  }
528  }
529  }
530  }
531  }
532 
533  if (cfg.StoreSDParticles()) {
534  const auto& pTypes = cfg.StoreSDParticleTypes();
535  for (const auto& p : simStat.ParticlesRange()) {
536  if (!pTypes.empty() && Is(std::abs(p.GetType())).In(pTypes))
537  continue;
538  SDParticle myParticle;
539  myParticle.SetType(p.GetType());
540  myParticle.SetEnergy(p.GetTotalEnergy());
541  myParticle.SetMomentum(ToTVector3(p.GetMomentum(), referenceCS /* what's adst unit for momentum? */));
542  myParticle.SetPosition(ToTVector3(p.GetPosition(), referenceCS, meter));
543  myParticle.SetTime(p.GetTime().GetSecond(), p.GetTime().GetNanoSecond());
544  myStation.AddParticle(myParticle);
545  }
546  }
547 
548  adstSEvent.AddSimStation(myStation);
549 
550  }
551 }
552 
553 
554 void
555 SD2ADST::FillUniversality(const evt::Event& event, RecEvent& recEvent)
556  const
557 {
558  auto& adstSEvent = recEvent.GetSDEvent();
559  auto& adstUnivRecShower= adstSEvent.GetUnivRecShower();
560 
561  if (!event.HasRecShower() || !event.GetRecShower().HasUnivRecShower())
562  return;
563 
564  const auto& showerUnivRecData = event.GetRecShower().GetUnivRecShower();
565  const auto& referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
566 
567  const auto& core = showerUnivRecData.GetCorePosition();
568  const auto& localCS = LocalCoordinateSystem::Create(core);
569 
570  double northing = 0;
571  double easting = 0;
572  double altitude = 0;
573 
574  try {
575  const utl::UTMPoint UTM(core, ReferenceEllipsoid::GetWGS84());
576  northing = UTM.GetNorthing();
577  easting = UTM.GetEasting();
578  altitude = UTM.GetHeight();
579  } catch (UTMPoint::UTMZoneException&) {
580  ERROR ("UTMZoneException: univ rec shower invalid");
581  return;
582  } catch (UTMPoint::UTMException&) {
583  ERROR ("UTMException: univ rec shower invalid");
584  return;
585  }
586 
587  // should be extended in the future to different rec types and so on...
588  adstUnivRecShower.SetRecLevel(showerUnivRecData.GetGoodRec() ? eGoodReconstruction: eNoReconstruction);
589 
590  adstUnivRecShower.SetCoreSiteCS(ToTVector3(core, referenceCS, meter));
591  adstUnivRecShower.SetCoreUTMCS(TVector3(easting, northing, altitude));
592 
593  const auto& coreError = showerUnivRecData.GetCoreError();
594  adstUnivRecShower.SetCoreNorthingError(coreError.GetY(localCS));
595  adstUnivRecShower.SetCoreEastingError(coreError.GetX(localCS));
596 
597  adstUnivRecShower.SetCoreTime(showerUnivRecData.GetCoreTime().GetGPSSecond(),
598  showerUnivRecData.GetCoreTime().GetGPSNanoSecond());
599 
600  const auto& axis = showerUnivRecData.GetAxis();
601 
602  TVector3 axisCore(0,0,1);
603  axisCore.SetTheta(axis.GetTheta(localCS));
604  axisCore.SetPhi(axis.GetPhi(localCS));
605  adstUnivRecShower.SetAxisCoreCS(axisCore);
606  adstUnivRecShower.SetZenithError(showerUnivRecData.GetThetaError());
607  adstUnivRecShower.SetAzimuthError(showerUnivRecData.GetPhiError());
608 
609  TVector3 axisSite(0,0,1);
610  axisSite.SetTheta(axis.GetTheta(referenceCS));
611  axisSite.SetPhi(axis.GetPhi(referenceCS));
612  adstUnivRecShower.SetAxisSiteCS(axisSite);
613 
614  adstUnivRecShower.SetEnergy(showerUnivRecData.GetEnergy());
615  adstUnivRecShower.SetEnergyError(showerUnivRecData.GetEnergyError());
616 
617  adstUnivRecShower.SetNmu(showerUnivRecData.GetNmu());
618  adstUnivRecShower.SetNmuError(showerUnivRecData.GetNmuError());
619 
620  const double gcm2 = g/cm2;
621 
622  adstUnivRecShower.SetXmax(showerUnivRecData.GetXmax()/gcm2);
623  adstUnivRecShower.SetXmaxError(showerUnivRecData.GetXmaxError()/gcm2);
624 
625  adstUnivRecShower.SetXmaxMu(showerUnivRecData.GetXmaxMu()/gcm2);
626  adstUnivRecShower.SetXmaxMuError(showerUnivRecData.GetXmaxMuError()/gcm2);
627 
628  adstUnivRecShower.SetLnA(showerUnivRecData.GetLnA());
629  adstUnivRecShower.SetLnAError(showerUnivRecData.GetLnAError());
630 
631  adstUnivRecShower.SetTimeModelOffset(showerUnivRecData.GetTimeModelOffset());
632  adstUnivRecShower.SetTimeModelOffsetError(showerUnivRecData.GetTimeModelOffsetError());
633 
634  adstUnivRecShower.SetNumberOfShapeCandidates(showerUnivRecData.GetNumberOfShapeCandidates());
635  adstUnivRecShower.SetNumberOfStartTimeCandidates(showerUnivRecData.GetNumberOfStartTimeCandidates());
636  adstUnivRecShower.SetNumberOfLDFCandidates(showerUnivRecData.GetNumberOfLDFCandidates());
637 
638  FillCelestialCoordinates(adstUnivRecShower);
639 }
640 
641 
642 void
643 SD2ADST::FillSEvent(const evt::Event& event, RecEvent& recEvent)
644  const
645 {
646  auto& adstSEvent = recEvent.GetSDEvent();
647  const auto& sevent = event.GetSEvent();
648  const auto& sDetector = det::Detector::GetInstance().GetSDetector();
649  auto& adstSdRecShower = adstSEvent.GetSdRecShower();
650  adstSEvent.SetRecLevel(eNoSdEvent);
651 
652  // the sd header info
653  adstSEvent.SetId(sevent.GetHeader().GetId());
654  const auto& eventTime = sevent.GetHeader().GetTime();
655  adstSEvent.SetEventTime(eventTime.GetGPSSecond(), eventTime.GetGPSNanoSecond());
656  const int yymmdd = otoa::TimeStamp2YYMMDD(eventTime);
657  const int hhmmss = otoa::TimeStamp2HHMMSS(eventTime);
658  adstSEvent.SetYYMMDD(yymmdd);
659  adstSEvent.SetHHMMSS(hhmmss);
660 
661  int nCandidates = 0;
662  int nAccidentals = 0;
663  double gpsSum = 0;
664  for (const auto& s : sevent.StationsRange()) {
665  if (s.IsRejected())
666  ++nAccidentals;
667  else if (s.IsCandidate()) {
668  ++nCandidates;
669  gpsSum += sDetector.GetStation(s.GetId()).GetCommissionTime().GetGPSSecond();
670  }
671  }
672 
673  adstSEvent.SetNumberOfAccidentalStations(nAccidentals);
674  adstSEvent.SetNumberOfCandidates(nCandidates);
675  adstSEvent.SetAverageStationAge(eventTime.GetGPSSecond() - gpsSum/nCandidates);
676 
677  if (nCandidates != 0)
678  adstSEvent.SetRecLevel(eHasTriggeredStations);
679 
680  if (sevent.HasTrigger()) {
681  adstSEvent.SetCDASSender(sevent.GetTrigger().GetSender());
682  adstSEvent.SetCDASAlgorithm(sevent.GetTrigger().GetAlgorithm());
683  }
684 
685  if (!event.HasRecShower() || !event.GetRecShower().HasSRecShower())
686  return;
687 
688  const auto& showerSRecData = event.GetRecShower().GetSRecShower();
689 
690  adstSEvent.SetBadPeriodId(showerSRecData.GetBadPeriodId());
691 
692  if (showerSRecData.IsT4())
693  adstSEvent.SetT4Trigger(showerSRecData.GetT4Trigger());
694  if (showerSRecData.IsT5())
695  adstSEvent.SetT5Trigger(showerSRecData.GetT5Trigger());
696  adstSEvent.SetT5HottestStation(showerSRecData.GetT5HottestStation());
697  adstSEvent.SetT5ClosestStation(showerSRecData.GetT5ClosestStation());
698  adstSEvent.SetT5ClosestStationIsUUB(showerSRecData.IsT5ClosestStationUUB());
699  adstSEvent.SetT5PriorActiveNeighbors(showerSRecData.GetT5PriorActiveNeighbors());
700  adstSEvent.SetT5PriorActiveUUBNeighbors(showerSRecData.GetT5PriorActiveUUBNeighbors());
701  adstSEvent.SetT5PostActiveNeighbors(showerSRecData.GetT5PosteriorActiveNeighbors());
702  adstSEvent.SetT5PostActiveUUBNeighbors(showerSRecData.GetT5PosteriorActiveUUBNeighbors());
703  adstSEvent.SetT5PostCoreTriangle(showerSRecData.GetT5PosteriorCoreTriangle());
704 
705  const auto& referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
706 
707  // seed
708  adstSdRecShower.SetSeedStations(showerSRecData.GetReconstructionSeed());
709  adstSdRecShower.SetSeedAxis(ToTVector3(showerSRecData.GetSeedAxis(), referenceCS));
710 
711  // core
712  const auto& core = showerSRecData.GetCorePosition();
713  const auto& localCS = LocalCoordinateSystem::Create(core);
714 
715  double northing = 0;
716  double easting = 0;
717  double altitude = 0;
718 
719  try {
720  utl::UTMPoint UTM(core, ReferenceEllipsoid::GetWGS84());
721  northing = UTM.GetNorthing();
722  easting = UTM.GetEasting();
723  altitude = UTM.GetHeight();
724  } catch (UTMPoint::UTMZoneException&) {
725  ERROR ("UTMZoneException: sd rec shower invalid");
726  return;
727  } catch (UTMPoint::UTMException&) {
728  ERROR ("UTMException: sd rec shower invalid");
729  return;
730  }
731 
732  adstSdRecShower.SetCoreSiteCS(ToTVector3(core, referenceCS, meter));
733  adstSdRecShower.SetCoreUTMCS(TVector3(easting, northing, altitude));
734 
735  const auto& coreError = showerSRecData.GetCoreError();
736  adstSdRecShower.SetCoreNorthingError(coreError.GetY(localCS));
737  adstSdRecShower.SetCoreEastingError(coreError.GetX(localCS));
738  adstSdRecShower.SetCoreNorthingEastingCorrelation(showerSRecData.GetCorrelationXY());
739 
740  adstSdRecShower.SetCoreTime(showerSRecData.GetCoreTime().GetGPSSecond(),
741  showerSRecData.GetCoreTime().GetGPSNanoSecond());
742 
743  // axis
744  adstSEvent.SetRecLevel(eHasSdAxis);
745  const auto& axis = showerSRecData.GetAxis();
746 
747  TVector3 axisCore(0,0,1);
748  axisCore.SetTheta(axis.GetTheta(localCS));
749  axisCore.SetPhi(axis.GetPhi(localCS));
750  adstSdRecShower.SetAxisCoreCS(axisCore);
751  adstSdRecShower.SetZenithError(showerSRecData.GetThetaError());
752  adstSdRecShower.SetAzimuthError(showerSRecData.GetPhiError());
753  adstSdRecShower.SetZenithAzimuthCorrelation(showerSRecData.GetCorrelationThetaPhi());
754 
755  TVector3 axisSite(0,0,1);
756  axisSite.SetTheta(axis.GetTheta(referenceCS));
757  axisSite.SetPhi(axis.GetPhi(referenceCS));
758  adstSdRecShower.SetAxisSiteCS(axisSite);
759 
760  adstSdRecShower.SetAngleChi2(showerSRecData.GetAngleChi2(), showerSRecData.GetAngleNdof());
761  adstSdRecShower.SetTimeResidualMean(showerSRecData.GetTimeResidualMean());
762  adstSdRecShower.SetTimeResidualSpread(showerSRecData.GetTimeResidualSpread());
763 
764  // save data of the true plane front reconstruction
765  if (showerSRecData.HasPlaneFrontRecData()) {
766  adstSdRecShower.SetPlaneFrontAxis(
767  ToTVector3(showerSRecData.GetPlaneFrontRecData().GetAxis(), referenceCS)
768  );
769  }
770 
771  // LDF and energy
772  if (showerSRecData.HasLDF()) {
773 
774  adstSEvent.SetRecLevel(eHasLDF);
775 
776  auto& adstLDF = adstSdRecShower.GetLDF();
777 
778  const double showerSize = showerSRecData.GetShowerSize();
779  const auto& ldfTable = showerSRecData.GetLDF();
780 
781  const size_t n = ldfTable.GetNPoints();
782  vector<double> rhos;
783  rhos.reserve(n);
784  vector<double> values;
785  values.reserve(n);
786  for (size_t i = 0; i < n; ++i) {
787  const double x = ldfTable[i].X();
788  rhos.push_back(x);
789  const double y = ldfTable[i].Y() / showerSize;
790  values.push_back(y);
791  }
792  adstLDF.SetTabulatedValues(rhos, values);
793 
794  adstLDF.SetLDFChi2(showerSRecData.GetLDFChi2(), showerSRecData.GetLDFNdof());
795  adstLDF.SetLDFLikelihood(abs(showerSRecData.GetLDFLikelihood()));
796 
797  adstLDF.SetBeta(showerSRecData.GetBeta(), showerSRecData.GetBetaError());
798  adstLDF.SetBetaSystematics(showerSRecData.GetBetaSystematics());
799  adstLDF.SetGamma(showerSRecData.GetGamma(), showerSRecData.GetGammaError());
800  adstLDF.SetNKGFermiParameters(showerSRecData.GetNKGFermiMu(),
801  showerSRecData.GetNKGFermiTau());
802 
803  adstLDF.SetLDFStatus(showerSRecData.GetLDFRecStage());
804 
805  const char* const showerSizeLabel = showerSRecData.GetShowerSizeLabel();
806 
807  if (showerSizeLabel && showerSizeLabel[0] == 'S') {
808  const unsigned int rRef = boost::lexical_cast<unsigned int>(showerSizeLabel + 1);
809  adstLDF.SetReferenceDistance(rRef);
810  } else
811  adstLDF.SetReferenceDistance(0);
812 
813  adstLDF.SetShowerSizeLabel(showerSizeLabel);
814  adstLDF.SetShowerSize(showerSize, showerSRecData.GetShowerSizeError());
815  adstLDF.SetShowerSizeSystematics(showerSRecData.GetShowerSizeSystematics());
816  adstLDF.SetSizeGeomagneticCorr(showerSRecData.GetSizeGeomagneticCorr(),
817  showerSRecData.GetSizeGeomagneticCorrError());
818  adstLDF.SetSizeWeatherCorr(showerSRecData.GetSizeAtmosphericCorr(),
819  showerSRecData.GetSizeAtmosphericCorrError());
820  }
821 
822  // Adding shower variables from scintillator
823  if (showerSRecData.HasScintillatorRecShower()) {
824 
825  const auto& showerScintillatorRecData = showerSRecData.GetScintillatorRecShower();
826  auto& adstScintillatorRecShower = adstSdRecShower.GetScintillatorRecShower();
827 
828  if (showerScintillatorRecData.HasLDF()) {
829 
830  auto& adstScintillatorLDF = adstScintillatorRecShower.GetLDF();
831  const double showerSize = showerScintillatorRecData.GetShowerSize();
832 
833  adstScintillatorLDF.SetLDFChi2(showerScintillatorRecData.GetLDFChi2(), showerScintillatorRecData.GetLDFNdof());
834  adstScintillatorLDF.SetLDFLikelihood(abs(showerScintillatorRecData.GetLDFLikelihood()));
835  adstScintillatorLDF.SetBeta(showerScintillatorRecData.GetBeta(), showerScintillatorRecData.GetBetaError());
836  adstScintillatorLDF.SetGamma(showerScintillatorRecData.GetGamma(), showerScintillatorRecData.GetGammaError());
837  adstScintillatorLDF.SetShowerSize(showerSize, showerScintillatorRecData.GetShowerSizeError());
838  adstScintillatorLDF.SetLDFStatus(showerScintillatorRecData.GetLDFRecStage());
839 
840  }
841 
842  }
843 
844  adstSdRecShower.SetROpt(showerSRecData.GetLDFOptimumDistance(), 0);
845 
846  adstSdRecShower.SetCicRefAngle(showerSRecData.GetCicRefAngle());
847  adstSdRecShower.SetAttShowerSize(showerSRecData.GetS38());
848 
849  const double energy = showerSRecData.GetEnergy();
850  adstSdRecShower.SetEnergy(energy);
851  adstSdRecShower.SetEnergyError(showerSRecData.GetEnergyError());
852  adstSdRecShower.SetEnergyLDFSys(showerSRecData.GetEnergyLDFSystematics());
853  adstSdRecShower.SetCurvature(showerSRecData.GetCurvature(), showerSRecData.GetCurvatureError());
854  adstSdRecShower.SetCurvatureOffset(showerSRecData.GetCurvatureTimeOffset(), showerSRecData.GetCurvatureTimeOffsetError());
855  adstSdRecShower.SetRadiusPoint(ToTVector3(showerSRecData.GetShowerCenter(), referenceCS, meter));
856 
857  // rise time related variables
858  if (showerSRecData.HasRiseTime1000()) {
859  const auto& sdrise = showerSRecData.GetRiseTime1000();
860  adstSdRecShower.GetRiseTimeResults().SetRiseTime1000(sdrise.GetRiseTime1000(),
861  sdrise.GetRiseTime1000Error());
862  adstSdRecShower.GetRiseTimeResults().SetXmax(sdrise.GetXmax());
863  adstSdRecShower.GetRiseTimeResults().SetXmaxErrorUp(sdrise.GetXmaxErrorUp());
864  adstSdRecShower.GetRiseTimeResults().SetXmaxErrorDown(sdrise.GetXmaxErrorDown());
865  adstSdRecShower.GetRiseTimeResults().SetRiseTimeChi2(sdrise.GetChi2());
866  adstSdRecShower.GetRiseTimeResults().SetRiseTimeNDF(sdrise.GetNdof());
867  adstSdRecShower.GetRiseTimeResults().SetAlpha(sdrise.GetAlpha());
868  adstSdRecShower.GetRiseTimeResults().SetBeta(sdrise.GetBeta());
869  }
870 
871  if (showerSRecData.HasFootprintData()) {
872  const auto& footprint = showerSRecData.GetFootprintData();
873  SdFootprintData& adstFootprint = adstSdRecShower.GetFootprint();
874  adstFootprint.SetLength(footprint.GetLength());
875  adstFootprint.SetWidth(footprint.GetWidth());
876  adstFootprint.SetAreaOverPeakAsymmetry(footprint.GetAreaOverPeakAsymmetry());
877  adstFootprint.SetSpeed(footprint.GetSpeed(), footprint.GetSpeedStandardDeviation());
878  adstFootprint.SetTOTFraction(footprint.GetTOTFraction());
879 
880  const auto& params = footprint.GetParameterVector();
881  for (unsigned int i = 0, n = params.size(); i < n; ++i)
882  adstFootprint.SetParameter(params[i], footprint.GetParameter(params[i]));
883  }
884 
885  FillCelestialCoordinates(adstSdRecShower);
886 
887  if (showerSRecData.GetCurvature())
888  adstSEvent.SetRecLevel(eHasCurvature);
889 
890  // filling meteo weather data
891  if (sevent.HasMeteo()) {
892  const auto& meteo = sevent.GetMeteo();
893  const auto& press = meteo.GetPressures();
894  const auto& temps = meteo.GetTemperatures();
895  const auto& humids = meteo.GetHumidities();
896  const auto& daypress = meteo.GetDayPressures();
897  const auto& daytemps = meteo.GetDayTemperatures();
898  const auto& dayhumids = meteo.GetDayHumidities();
899 
900  // assuming that the fields have the same length, should be safe
901  for (unsigned int i = 0, n = press.size(); i < n; ++i) {
902 
903  MeteoData idata;
904  idata.SetWeatherInformation(ePressure, press[i]/hPa);
905  idata.SetWeatherInformation(eTemperature, temps[i]);
906  idata.SetWeatherInformation(eHumidity, humids[i]);
907  idata.SetWeatherInformation(eDayPressure, daypress[i]/hPa);
908  idata.SetWeatherInformation(eDayTemperature, daytemps[i]);
909  idata.SetWeatherInformation(eDayHumidity, dayhumids[i]);
910 
911  recEvent.GetDetector().AddMeteoSiteInfo(idata);
912  }
913  }
914  //Parameters from Parameter Storage
915  const auto& sRecShower = event.GetRecShower().GetSRecShower();
916  for (const auto par : sRecShower.GetEnumVector()) {
917  if (sRecShower.HasParameter(par))
918  adstSdRecShower.SetParameter(par, sRecShower.GetParameter(par));
919  else
920  std::cout << "WARNING: Parameter " << par << " does not exist" << std::endl;
921  }
922 
923  int ncount = 0;
924  for (const auto& ij : sRecShower.GetCovarianceEnumVector()) {
925  const auto& i = ij.first;
926  const auto& j = ij.second;
927  if (sRecShower.HasParameterCovariance(i, j))
928  adstSdRecShower.SetParameterCovariance(i, j, sRecShower.GetParameterCovariance(i, j));
929  else
930  ++ncount;
931  }
932 
933  // for debugging purpose, will be removed later
934  if (ncount)
935  std::cout << "WARNING: " << ncount << " ParameterCovariance pairs does not exist" << endl;
936 }
937 
938 
939 void
940 SD2ADST::FillMPD(const evt::Event& event, RecEvent& recEvent)
941  const
942 {
943  if (event.HasRecShower() &&
944  event.GetRecShower().HasSRecShower()) {
945  const auto& sdRecShower = event.GetRecShower().GetSRecShower();
946  if (sdRecShower.HasMPDHistogram()) {
947  const auto& mpd = sdRecShower.GetMPDHistogram();
948  auto& sd = recEvent.GetSDEvent();
949  sd.GetSdRecShower().SetMPDMin(mpd.GetStart());
950  sd.GetSdRecShower().SetMPDMax(mpd.GetStop());
951  sd.GetSdRecShower().SetMPDData(mpd.GetBins());
952  }
953  }
954 }
955 
956 
957 ETraceType
958 SD2ADST::GetTraceType(const sevt::StationConstants::SignalComponent comp)
959  const
960 {
961  switch (comp) {
962  case sevt::StationConstants::eTotal: return eTotalTrace;
963  case sevt::StationConstants::eTotalNoSaturation: return eTotalUnsaturatedTrace;
964  case sevt::StationConstants::ePhoton: return ePhotonTrace;
965  case sevt::StationConstants::eElectron: return eElectronTrace;
966  case sevt::StationConstants::eMuon: return eMuonTrace;
967  case sevt::StationConstants::eShowerMuonDecayPhoton: return ePhotonFromMuonTrace;
968  case sevt::StationConstants::eShowerMuonDecayElectron: return eElectronFromMuonTrace;
969  case sevt::StationConstants::eHadron: return eHadronTrace;
970  // components for universality parametrization
971  case sevt::StationConstants::eShowerLocalHadronPhoton: return eShowerLocalHadronPhotonTrace;
972  case sevt::StationConstants::eShowerLocalHadronElectron: return eShowerLocalHadronElectronTrace;
973  default: return eUnknownTrace;
974  }
975 }
976 
977 
978 #ifdef INCLUDE_HAS
979 void
980 RecDataWriter::FillMuonMaps(const RecDataProvider& theRec, SDEvent& theSd)
981 {
982  if (theRec.GetSdLDFStatus() < 13) {
983  INFO("No SdHorizontalReconstruction. No muon map written.");
984  return;
985  }
986 
987  const auto& referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
988 
989  if (!fMuonProfile || !fTankResponse || !fEMComponent) {
990  INFO("SdHorizontalReconstruction done, but RecDataWriter not fully configured. No muon map written.");
991  return;
992  }
993 
994  const double theta = theRec.GetSdTheta();
995  const double phi = theRec.GetSdPhi();
996  const double n19 = theRec.GetSdS1000();
997  const double xCore = theRec.GetSdCoreSite().get<0>();
998  const double yCore = theRec.GetSdCoreSite().get<1>();
999  const int squarerootRBinsInScan = 200; //63;
1000  const double outerRadiusMax = squarerootRBinsInScan*squarerootRBinsInScan;
1001 
1002  MuonMap& myMap = theSd.GetMuonMap();
1003 
1004  unsigned int sqrtRStop[fNoOfPoints];
1005  for (unsigned int phiIter = 0; phiIter < fNoOfPoints; ++phiIter)
1006  sqrtRStop[phiIter] = 1; // remember radius of last contour
1007 
1008  for (int contourIter = fNoOfContours; contourIter > 0; --contourIter) {
1009  const double myContourSignal = fContourSignalFactor * sqrt(n19) * exp(0.5 * contourIter);
1010 
1011  MuonMapContour myContour;
1012  myContour.SetContourSignal(myContourSignal);
1013 
1014  unsigned int nrOfPointsNotSet = 0;
1015 
1016  for (unsigned int phiIter = 0; phiIter < fNoOfPoints; ++phiIter) {
1017  const double localPhi = phiIter * 2 * kPi / fNoOfPoints;
1018  double outerRadius = 0;
1019  double innerRadius = 1;
1020  double contourRadius = 0;
1021 
1022  double outerMuonSignal = 0;
1023  double innerMuonSignal =
1024  GetMuonMapSignal(innerRadius * cos(localPhi), innerRadius * sin(localPhi), theta, phi, n19);
1025 
1026  for (int sqrtRIter = sqrtRStop[phiIter]; sqrtRIter < squarerootRBinsInScan; ++sqrtRIter) {
1027  innerRadius = outerRadius;
1028  outerRadius = sqrtRIter*sqrtRIter;
1029  innerMuonSignal = outerMuonSignal;
1030  outerMuonSignal =
1031  GetMuonMapSignal(outerRadius * cos(localPhi), outerRadius * sin(localPhi), theta, phi, n19);
1032  if (outerMuonSignal < 0.2)
1033  contourRadius = outerRadius;
1034  else if (myContourSignal > innerMuonSignal || outerMuonSignal > myContourSignal)
1035  continue;
1036  else if (innerMuonSignal >= myContourSignal && myContourSignal >= outerMuonSignal)
1037  contourRadius = // interpolate!
1038  ((innerRadius - outerRadius) / (innerMuonSignal - outerMuonSignal)) * (myContourSignal - outerMuonSignal) + outerRadius;
1039  sqrtRStop[phiIter] = sqrtRIter - 1;
1040  break;
1041  }
1042  if (contourRadius > 0) {
1043  const double electronSignal = myContourSignal * fEMComponent->SignalRatio(contourRadius * cos(localPhi),
1044  contourRadius * sin(localPhi), theta, phi);
1045  const double xGroundSite = xCore + contourRadius * cos(localPhi);
1046  const double yGroundSite = yCore + contourRadius * sin(localPhi);
1047  const MuonMapPoint myPoint(xGroundSite, yGroundSite, electronSignal);
1048  myContour.AddPoint(myPoint);
1049  if (fVerbosity > 1) {
1050  ostringstream msg;
1051  msg << "For muon contour at " << myContourSignal << " VEM: Found point at r = " << contourRadius << ".";
1052  INFO(msg);
1053  }
1054  } else {
1055  ++nrOfPointsNotSet;
1056  if (fVerbosity > 1) {
1057  ostringstream msg;
1058  msg << "For muon contour at " << myContourSignal << " VEM: " << nrOfPointsNotSet << " points of contour not found within range.";
1059  INFO(msg);
1060  }
1061  }
1062  }
1063  if (nrOfPointsNotSet <= (unsigned int)(0.5 * fNoOfPoints))
1064  myMap.AddContour(myContour);
1065  else {
1066  if (fVerbosity > 0) {
1067  ostringstream msg;
1068  msg << "Contour at " << myContourSignal << " VEM not added to map: Not enough points (" << fNoOfPoints-nrOfPointsNotSet
1069  << " of " << fNoOfPoints << ") in range.";
1070  INFO(msg);
1071  }
1072  }
1073  }
1074 
1075  myMap.SetMuonMapType(fMuonMapType);
1076  myMap.SetMagneticFieldAzimuth(HASUtilities::GetBpsi(theta, phi));
1077 }
1078 
1079 
1080 double
1081 RecDataWriter::GetMuonMapSignal(double x, double y, double theta, double phi, double n19)
1082 {
1083  double nMu = 0;
1084  try {
1085  nMu = n19*fMuonProfile->NMuon(x, y, theta, phi, 10*EeV);
1086  } catch (OutOfBoundException&) {
1087  return 0;
1088  }
1089  if (isinf(nMu) || std::isnan(nMu)) {
1090  INFO("Unphysical predicted muon number!");
1091  return 0;
1092  }
1093  int nMuDown = int(nMu);
1094  int nMuUp = int(nMu) + 1;
1095 
1096  const double precision = 1e-2;
1097  double sum = 0;
1098 
1099  for (int i = 0; i < 100; ++i) {
1100  double val = 0;
1101  if (nMuDown > 0) {
1102  val = TMath::Poisson(nMuDown, nMu) * fTankResponse->Moment(theta, nMuDown, 1);
1103  --nMuDown;
1104  }
1105  val += TMath::Poisson(nMuUp, nMu) * fTankResponse->Moment(theta, nMuUp, 1);
1106  ++nMuUp;
1107  if (val < precision*sum || val == 0)
1108  break;
1109  sum += val;
1110  }
1111 
1112  return sum;
1113 }
1114 
1115 
1116 double
1117 RecDataWriter::GetMuonMapError(double x, double y, double theta, double phi, double n19)
1118 {
1119  double nMu = 0;
1120  try {
1121  nMu = n19*fMuonProfile->NMuon(x, y, theta, phi, 10*EeV);
1122  } catch (OutOfBoundException& e) {
1123  return 0;
1124  }
1125  if (isinf(nMu) || std::isnan(nMu)) {
1126  INFO("Unphysical predicted muon number!");
1127  return 0;
1128  }
1129  int nMuDown = int(nMu);
1130  int nMuUp = int(nMu) + 1;
1131 
1132  const double precision = 1e-2;
1133  double sum = 0;
1134 
1135  for (int i = 0; i < 100; ++i) {
1136  double val = 0;
1137  if (nMuDown > 0) {
1138  const double p = TMath::Poisson(nMuDown, nMu);
1139  const double mom1 = fTankResponse->Moment(theta, nMuDown, 1);
1140  const double mom2 = fTankResponse->Moment(theta, nMuDown, 2);
1141  val = p * sqrt(mom2 - mom1*mom1);
1142  --nMuDown;
1143  }
1144  const double p = TMath::Poisson(nMuUp, nMu);
1145  const double mom1 = fTankResponse->Moment(theta, nMuUp, 1);
1146  const double mom2 = fTankResponse->Moment(theta, nMuUp, 2);
1147  val += p * sqrt(mom2 - mom1*mom1);
1148  ++nMuUp;
1149 
1150  if (val < precision*sum || val == 0)
1151  break;
1152 
1153  sum += val;
1154  }
1155 
1156  sum *= 1 + fEMComponent->SignalRatio(x, y, theta, phi);
1157 
1158  return sum;
1159 }
1160 
1161 
1162 #endif
io::Convert
void Convert(std::vector< T1, A1 > &destination, const std::vector< T2, A2 > &source)
Definition: StreamerUtilities.h:32
utl::Is
IsProxy< T > Is(const T &x)
Definition: Is.h:46
sevt::StationConstants::eTotal
total (shower and background)
Definition: SEvent/StationConstants.h:13
utl::AugerException
Base class for all exceptions used in the auger offline code.
Definition: AugerException.h:24
evt::Event::HasRecShower
bool HasRecShower() const
Definition: Framework/Event/Event.h:55
utl::UTMPoint
Class to hold and convert a point in geodetic coordinates.
Definition: UTMPoint.h:40
evt::Event::HasSimShower
bool HasSimShower() const
Definition: Framework/Event/Event.h:50
tls::GetBpsi
double GetBpsi(const double theta, const double phi)
Definition: MuonProfileUtilities.h:19
galactic::meter
const double meter
Definition: GalacticUnits.h:29
INFO
#define INFO(message)
Macro for logging informational messages.
Definition: ErrorLogger.h:161
sevt::StationConstants::eShowerLocalHadronElectron
electrons produced by hadrons close to the detector
Definition: SEvent/StationConstants.h:65
otoa::TimeStamp2HHMMSS
unsigned int TimeStamp2HHMMSS(const utl::TimeStamp &timest)
Convert a TimeStamp into an integer representing the time as HHMMSS.
Definition: ConversionUtil.cc:46
gain
int gain
Definition: dump1090.h:241
galactic::EeV
const double EeV
Definition: GalacticUnits.h:34
utl::UTMPoint::GetNorthing
double GetNorthing() const
Get the northing.
Definition: UTMPoint.h:206
utl::UTMPoint::UTMZoneException
Report attempts to use invalid UTM zone.
Definition: UTMPoint.h:300
utl::TimeStamp
A TimeStamp holds GPS second and nanosecond for some event.
Definition: TimeStamp.h:110
utl::OutOfBoundException
Exception for reporting variable out of valid range.
Definition: AugerException.h:118
sevt::StationConstants::eElectron
electrons and positrons from shower
Definition: SEvent/StationConstants.h:15
sevt::ByDecreasingSignal
Criterion to sort stations by decreasing signal.
Definition: SortCriteria.h:41
precision
static const double precision
Definition: VTankResponse.cc:12
sevt::StationConstants::eTotalNoSaturation
total FADC trace, with no saturation applied by FADC/baseline simulator
Definition: SEvent/StationConstants.h:53
utl::s
constexpr double s
Definition: AugerUnits.h:163
utl::UTMPoint::GetHeight
double GetHeight() const
Get the height.
Definition: UTMPoint.h:212
utl::kPi
constexpr double kPi
Definition: MathConstants.h:24
utl::nanosecond
constexpr double nanosecond
Definition: AugerUnits.h:143
utl::abs
double abs(const SVector< n, T > &v)
Definition: ValidationTests/SdCalibration/Validatrix.cc:32
evt::Event::GetSimShower
ShowerSimData & GetSimShower()
Definition: Framework/Event/Event.cc:31
utl::UTMPoint::GetEasting
double GetEasting() const
Get the easting.
Definition: UTMPoint.h:209
utl::g
constexpr double g
Definition: AugerUnits.h:200
otoa::TimeStamp2YYMMDD
unsigned int TimeStamp2YYMMDD(const utl::TimeStamp &timest)
Convert a TimeStamp into an integer representing the date as YYMMDD.
Definition: ConversionUtil.cc:38
utl::hPa
constexpr double hPa
Definition: AugerUnits.h:218
utl::kSpeedOfLight
constexpr double kSpeedOfLight
Definition: PhysicalConstants.h:59
utl::TimeInterval
A TimeInterval is used to represent time elapsed between two events.
Definition: TimeInterval.h:43
utl::UTMPoint::UTMException
Report problems in UTM handling.
Definition: UTMPoint.h:288
utl::Vector
Vector object.
Definition: Vector.h:30
gcm2
const double gcm2
Definition: LongitudinalXmaxScanner.cc:54
otoa::ToTVector3
TVector3 ToTVector3(const T &v, const utl::CoordinateSystemPtr &cs, const double unit=1)
Definition: ConversionUtil.h:51
otoa::FillCelestialCoordinates
void FillCelestialCoordinates(RecShower &recShower)
Definition: ConversionUtil.cc:181
ERROR
#define ERROR(message)
Macro for logging error messages.
Definition: ErrorLogger.h:165
utl::m
constexpr double m
Definition: AugerUnits.h:121
sevt::StationConstants::eMuon
mu+ and mu- (including signal from mu decay electrons) from shower
Definition: SEvent/StationConstants.h:19
utl::AugerException::GetMessage
const std::string & GetMessage() const
Retrieve the message from the exception.
Definition: AugerException.h:38
RPerp2
double RPerp2(const Vector &axis, const Vector &station)
Definition: SdSimMuonNumberFitter.cc:69
evt::ShowerSimData::GetShowerCoordinateSystem
utl::CoordinateSystemPtr GetShowerCoordinateSystem() const
Definition: ShowerSimData.h:227
GetConfig
std::string GetConfig()
Definition: fwkPython.cc:23
utl::cm2
constexpr double cm2
Definition: AugerUnits.h:118
, generated on Tue Sep 26 2023.