Modules/SdSimulation/Deprecated/UpgradeASCIITests/G4TankSimulatorASCII/DataWriter.cc
Go to the documentation of this file.
1 #include "DataWriter.h"
2 
3 #include <iostream>
4 #include <fstream>
5 #include <iterator>
6 
7 #include <det/Detector.h>
8 #include <sdet/SDetector.h>
9 #include <sdet/Station.h>
10 
11 #include <evt/Event.h>
12 
13 #include <fwk/RunController.h>
14 
15 #include <sevt/SEvent.h>
16 #include <sevt/Station.h>
17 #include <sevt/PMT.h>
18 #include <sevt/PMTRecData.h>
19 #include <sevt/PMTSimData.h>
20 #include <sevt/PMTCalibData.h>
21 #include <sevt/StationSimData.h>
22 #include <sevt/StationRecData.h>
23 
24 #include <utl/Trace.h>
25 #include <utl/AugerUnits.h>
26 #include <utl/MathConstants.h>
27 #include <utl/PhysicalConstants.h>
28 #include <utl/TimeStamp.h>
29 #include <utl/TabulatedFunction.h>
30 #include <utl/ErrorLogger.h>
31 #include <utl/Particle.h>
32 #include <utl/AugerCoordinateSystem.h>
33 #include <utl/MultiTimeDistribution.h>
34 
35 #include <TObjArray.h>
36 #include <TObject.h>
37 #include <TGraph.h>
38 #include <TGraphErrors.h>
39 #include <TMultiGraph.h>
40 #include <TMarker.h>
41 #include <TTree.h>
42 #include <TFolder.h>
43 #include <TString.h>
44 
45 #include <sstream>
46 #include <string>
47 #include <vector>
48 
49 #include "G4TankSimulator.h"
50 #include "G4TankTrackingAction.h"
51 
52 using namespace det;
53 using namespace evt;
54 using namespace fwk;
55 using namespace sevt;
56 using namespace utl;
57 using namespace G4TankSimulatorASCII;
58 
59 using namespace std;
60 
61 
62 VModule::ResultFlag DataWriter::Init(void) {
63  //INFO("DataWriter::Init()");
64 
65  PrintPulses=true;
66 
67  // Read File Name
68  std::string FileName="DataWriter.dat";
69 
70  fTextFile = new ofstream(FileName.c_str());
71  if ( !PrintPulses ) return eSuccess;
72  FileName="DataWriter.pulses";
73  fPulsesFile = new ofstream(FileName.c_str());
74 
75  return eSuccess;
76 }
77 
78 VModule::ResultFlag DataWriter::Run(Event &theEvent){
79 
80  //INFO("DataWriter::Run()");
81 
82  if (theEvent.HasSEvent())
83  WriteInfo(theEvent.GetSEvent());
84  else
85  INFO("DataWriter::No SEvent!!!!!!!!");
86  return eSuccess;
87 }
88 
89 VModule::ResultFlag DataWriter::Finish(void) {
90 
91  //INFO("DataWriter::Finish()");
92 
93  fTextFile->close();
94  if ( !PrintPulses ) return eSuccess;
95  fPulsesFile->close();
96 
97  return eSuccess;
98 
99 }
100 
101 
102 
103 void DataWriter::WriteInfo(sevt::SEvent &event){
104 
105  INFO("DataWriter::WriteInfo()");
106 
107  for (sevt::SEvent::StationIterator sIt = event.StationsBegin();
108  sIt != event.StationsEnd(); ++sIt) {
109 
110 
111  if (!sIt->HasSimData())
112  continue;
113 
114  sevt::StationSimData &simData = sIt->GetSimData();
115 
116  const unsigned long numParticles =
117  simData.ParticlesEnd() - simData.ParticlesBegin();
118 
119 
120  unsigned int fPESum[2]={0,0};// [ascin/wcd]
121  unsigned int fADCSum_HG[2]={0,0};
122  unsigned int fADCSum_LG[2]={0,0};
123 
124  //----
125  unsigned int npe=0;
126  for (sevt::Station::PMTIterator pmtIt = sIt->PMTsBegin(sdet::PMTConstants::eAnyType); pmtIt != sIt->PMTsEnd(sdet::PMTConstants::eAnyType); ++pmtIt )
127  {
128  if (!pmtIt->HasSimData()) continue;
129  if (!pmtIt->GetSimData().HasPETimeDistribution(sevt::StationConstants::eTotal)) continue;
130 
131  utl::TimeDistributionI &pe = pmtIt->GetSimData().GetPETimeDistribution(sevt::StationConstants::eTotal);
132  for (utl::TimeDistributionI::SparseIterator peIt =pe.SparseBegin();peIt != pe.SparseEnd();peIt++)
133  npe+=peIt->get<1>();
134  }
135 
136  if ( npe<=100 ) continue;
137 
138 
139  //----
140 
141  for (sevt::Station::PMTIterator pmtIt = sIt->PMTsBegin(sdet::PMTConstants::eAnyType); pmtIt != sIt->PMTsEnd(sdet::PMTConstants::eAnyType); ++pmtIt ) {
142 
143  if (pmtIt->GetType() == sdet::PMTConstants::eWaterCherenkovSmall) {
144  continue;
145  }
146 
147  int idetector = 1;
148  if (pmtIt->GetType() == sdet::PMTConstants::eScintillator) {
149  idetector = 0;
150  }
151  else if (pmtIt->GetType() != sdet::PMTConstants::eWaterCherenkovLarge) {
152  cout << "[ERROR]: Unrecognized PMT type. Terminating." << endl;
153  exit(1);
154  }
155 
156  sevt::StationConstants::SignalComponent component = sevt::StationConstants::eTotal;
157 
158  if (!pmtIt->HasSimData()) continue;
159  if (!pmtIt->GetSimData().HasPETimeDistribution(component)) continue;
160 
161 
162  //-- Photoelectron count
163  utl::TimeDistributionI &pe = pmtIt->GetSimData().GetPETimeDistribution(component);
164  double timeSlot= pe.GetBinning();
165  int entries=0;
166  for (utl::TimeDistributionI::SparseIterator peIt =pe.SparseBegin();peIt != pe.SparseEnd();peIt++)
167  {
168  entries++;
169  unsigned int npe=peIt->get<1>();
170  fPESum[idetector]+= npe;
171  }
172 
173  //-- ADC integral
174  if ( pmtIt->GetSimData().HasFADCTrace(component) )
175  {
176  for (int igain=0;igain<2;igain++)
177  {
178  TimeDistributionI& trace = ( igain==0 ? pmtIt->GetSimData().GetFADCTrace(PMTConstants::eHighGain,component) :
179  pmtIt->GetSimData().GetFADCTrace(PMTConstants::eLowGain, component) );
180 
181  for (int i = trace.GetStart(); i <= trace.GetStop(); ++i)
182  {
183  if ( igain==0 ) fADCSum_HG[idetector]+=trace[i];
184  else fADCSum_LG[idetector]+=trace[i];
185  }
186  }
187  }
188 
189  //-----
190  if ( !PrintPulses ) continue;
191 
192  //-- Print photoelectron trace
193 
194  if ( idetector == 0 ) *fPulsesFile << " PE_trace_ASCII " << sIt->GetId() << " " << entries << " " << endl;
195  if ( idetector == 1 && pmtIt->GetId()==1 ) *fPulsesFile << " PE_trace_WCD1 " << sIt->GetId() << " " << entries << " " << endl;
196  if ( idetector == 1 && pmtIt->GetId()==2 ) *fPulsesFile << " PE_trace_WCD2 " << sIt->GetId() << " " << entries << " " << endl;
197  if ( idetector == 1 && pmtIt->GetId()==3 ) *fPulsesFile << " PE_trace_WCD3 " << sIt->GetId() << " " << entries << " " << endl;
198 
199  if ( entries > 0 )
200  {
201  for (utl::TimeDistributionI::SparseIterator peIt =pe.SparseBegin();peIt != pe.SparseEnd();peIt++)
202  {
203  double time=timeSlot*(peIt->get<0>());
204  unsigned int npe=peIt->get<1>();
205  *fPulsesFile<< time <<" " << npe << " " ;
206  }
207  *fPulsesFile << endl;
208  }
209 
210  //-- Print FADC trace
211  if ( pmtIt->GetSimData().HasFADCTrace(component) )
212  {
213  for (int igain=0;igain<2;igain++)
214  {
215  TimeDistributionI& trace = ( igain==0 ? pmtIt->GetSimData().GetFADCTrace(PMTConstants::eHighGain,component) :
216  pmtIt->GetSimData().GetFADCTrace(PMTConstants::eLowGain, component) );
217 
218  int entries=0;
219  for (int i = trace.GetStart(); i <= trace.GetStop(); ++i)
220  if ( trace[i]>0 ) entries++;
221 
222  if ( idetector == 0 && igain==0 ) *fPulsesFile << " FADC_HG_ASCII " << sIt->GetId() << " " << entries << " " << endl;
223  if ( idetector == 0 && igain==1 ) *fPulsesFile << " FADC_LG_ASCII " << sIt->GetId() << " " << entries << " " << endl;
224  if ( idetector == 1 && igain==0 && pmtIt->GetId()==1 ) *fPulsesFile << " FADC_HG_WCD1 " << sIt->GetId() << " " << entries << " " << endl;
225  if ( idetector == 1 && igain==1 && pmtIt->GetId()==1 ) *fPulsesFile << " FADC_LG_WCD1 " << sIt->GetId() << " " << entries << " " << endl;
226  if ( idetector == 1 && igain==0 && pmtIt->GetId()==2 ) *fPulsesFile << " FADC_HG_WCD2 " << sIt->GetId() << " " << entries << " " << endl;
227  if ( idetector == 1 && igain==1 && pmtIt->GetId()==2 ) *fPulsesFile << " FADC_LG_WCD2 " << sIt->GetId() << " " << entries << " " << endl;
228  if ( idetector == 1 && igain==0 && pmtIt->GetId()==3 ) *fPulsesFile << " FADC_HG_WCD3 " << sIt->GetId() << " " << entries << " " << endl;
229  if ( idetector == 1 && igain==1 && pmtIt->GetId()==3 ) *fPulsesFile << " FADC_LG_WCD3 " << sIt->GetId() << " " << entries << " " << endl;
230 
231  if ( entries > 0 )
232  {
233  for (int i = trace.GetStart(); i <= trace.GetStop(); ++i)
234  if ( trace[i]>0 ) *fPulsesFile << i << " " << trace[i] << " " ;
235  *fPulsesFile << endl;
236  }
237  }
238  }
239  *fPulsesFile << endl;
240  }
241 
242 
243 
244  //---
245  //Output to generate the lookup tables
246  *fTextFile << sIt->GetId() << " " << fPESum[0] << " " << fPESum[1]
247  << " " << fADCSum_LG[0] << " " << fADCSum_HG[0] << " " << fADCSum_LG[1] << " " << fADCSum_HG[1]
248  << std::endl;
249 
250  ostringstream msg;
251  msg << "StatId= " << sIt->GetId()<< " " << fPESum[0] << " " << fPESum[1]
252  << " " << fADCSum_LG[0] << " " << fADCSum_HG[0] << " " << fADCSum_LG[1] << " " << fADCSum_HG[1]
253  << std::endl;
254  msg << "------------------------------------------------------" <<endl<< endl;
255  INFO(msg);
256 
257  }
258 
259 
260 }
261 
262 // Configure (x)emacs for this file ...
263 // Local Variables:
264 // modex:c++
265 // compile-command: "make -C .. -k"
266 // End:
DataWriter::Finish
fwk::VModule::ResultFlag Finish()
Finish: invoked at end of the run (NOT end of the event)
Definition: Modules/SdSimulation/Deprecated/UpgradeASCIITests/G4TankSimulatorASCII/DataWriter.cc:89
utl::TimeDistribution::GetStart
int GetStart() const
First slot with data.
Definition: TimeDistribution.h:207
sevt::StationSimData
Station Level Simulated Data
Definition: SEvent/StationSimData.h:32
sevt::StationConstants::eTotal
total (shower and background)
Definition: SEvent/StationConstants.h:13
sevt::SEvent
Interface class to access to the SD part of an event.
Definition: SEvent.h:39
utl::TimeDistribution
Histogram class for time distributions with suppressed empty bins.
Definition: TimeDistribution-fwd.h:15
INFO
#define INFO(message)
Macro for logging informational messages.
Definition: ErrorLogger.h:161
DataWriter::Init
fwk::VModule::ResultFlag Init()
Initialize: invoked at beginning of run (NOT beginning of event)
Definition: Modules/SdSimulation/Deprecated/UpgradeASCIITests/G4TankSimulatorASCII/DataWriter.cc:62
sevt::Station::PMTIterator
boost::filter_iterator< PMTFilter, InternalPMTIterator > PMTIterator
Iterator over station for read/write.
Definition: SEvent/Station.h:87
exit
int exit
Definition: dump1090.h:237
DataWriter::Run
fwk::VModule::ResultFlag Run(evt::Event &theEvent)
Run: invoked once per event.
Definition: Modules/SdSimulation/Deprecated/UpgradeASCIITests/G4TankSimulatorASCII/DataWriter.cc:78
utl::TimeDistribution::GetBinning
double GetBinning() const
Size of one slot.
Definition: TimeDistribution.h:204
utl::TimeDistribution::GetStop
int GetStop() const
Last slot with data (1 less than First slot if no data)
Definition: TimeDistribution.h:210
utl::TimeDistribution::SparseEnd
SparseIterator SparseEnd() const
Definition: TimeDistribution.h:219
sevt::StationSimData::ParticlesBegin
ParticleIterator ParticlesBegin()
Beginning of simulated particles entering the station.
Definition: SEvent/StationSimData.h:46
fwk::VModule::ResultFlag
ResultFlag
Flag returned by module methods to the RunController.
Definition: VModule.h:60
sevt::SEvent::StationsBegin
StationIterator StationsBegin()
Beginning of all stations.
Definition: SEvent.h:57
sevt::SEvent::StationIterator
boost::indirect_iterator< InternalStationIterator, Station & > StationIterator
Iterator over all stations.
Definition: SEvent.h:52
utl::TimeDistribution::SparseBegin
SparseIterator SparseBegin() const
Iterator over time slots with data in them (skips empty slots).
Definition: TimeDistribution.h:217
sevt::StationSimData::ParticlesEnd
ParticleIterator ParticlesEnd()
End of simulated particles entering the station.
Definition: SEvent/StationSimData.h:49
evt::Event::GetSEvent
sevt::SEvent & GetSEvent()
Definition: Framework/Event/Event.cc:81
evt::Event::HasSEvent
bool HasSEvent() const
Definition: Framework/Event/Event.h:65
utl::TimeDistribution::SparseIterator
boost::transform_iterator< InternalMapFunctor, InternalConstIterator, Tuple > SparseIterator
Definition: TimeDistribution.h:72
, generated on Tue Sep 26 2023.