source: trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc@ 9625

Last change on this file since 9625 was 9616, checked in by rohlfs, 14 years ago
can read also EventIO files
File size: 10.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19 Qi Zhe, 06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
20
21! Copyright: Software Development, 2000-2009
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MCorsikaRunHeader
29//
30// Root storage container for the RUN HEADER information
31//
32// Class Version 2:
33// ----------------
34// + UInt_t fParticleID
35// + Float_t fImpactMax
36// + Float_t fMagneticFieldX
37// + Float_t fMagneticFieldZ
38// + Float_t fMagneticFieldAz
39// + Float_t fAtmosphericLayers[5]
40// + Float_t fAtmosphericCoeffA[5]
41// + Float_t fAtmosphericCoeffB[5]
42// + Float_t fAtmosphericCoeffC[5]
43// + UInt_t fCerenkovFlag
44//
45////////////////////////////////////////////////////////////////////////////
46
47#include "MCorsikaRunHeader.h"
48#include "MCorsikaFormat.h"
49
50#include <fstream>
51#include <iomanip>
52
53#include "MLog.h"
54#include "MLogManip.h"
55
56#include "MMcEvt.hxx"
57
58ClassImp(MCorsikaRunHeader);
59
60using namespace std;
61
62const Double_t MCorsikaRunHeader::fgEarthRadius = 637131500; // [cm] Earth radius as defined in CORSIKA
63
64// --------------------------------------------------------------------------
65//
66// Default constructor. Creates array which stores the pixel assignment.
67//
68//
69MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title)
70 : fNumObsLevel(0), fImpactMax(-1), fZdMin(0), fZdMax(-1),
71 fAzMin(0), fAzMax(0), fViewConeInnerAngle(0), fViewConeOuterAngle(-1)
72{
73 fName = name ? name : "MCorsikaRunHeader";
74 fTitle = title ? title : "Raw Run Header Information";
75}
76
77// --------------------------------------------------------------------------
78//
79// Read in one run header from the binary file
80//
81Bool_t MCorsikaRunHeader::ReadEvt(MCorsikaFormat * fInFormat)
82{
83 if (!fInFormat->SeekNextBlock("RUNH", 1200))
84 return kFALSE;
85
86 Float_t f[272];
87 if (!fInFormat->ReadData(272, f))
88 return kFALSE;
89
90 fRunNumber = TMath::Nint(f[0]);
91 fNumEvents = 0;
92
93 fRunStart.SetCorsikaTime(f[1]);
94
95 fProgramVersion = f[2];
96 fNumObsLevel = TMath::Nint(f[3]);
97
98 if (fNumObsLevel!=1)
99 {
100 *fLog << err << "ERROR - Currently only one observation level is allowed." << endl;
101 return kFALSE;
102 }
103
104 memset(fObsLevel, 0, 10*4);
105 memcpy(fObsLevel, f+4, fNumObsLevel*4);
106
107 fSlopeSpectrum = f[14];
108 fEnergyMin = f[15];
109 fEnergyMax = f[16];
110
111 // Implemented in CORSIKA Version >= 6.822
112 fImpactMax = -1;
113
114 // CORSIKA scattering in a disc on the ground
115 if (f[246]>0 && f[247]==0)
116 {
117 *fLog << warn << "WARNING - Events scattered in a disc on the ground." << endl;
118 fImpactMax = f[246];
119 }
120
121 // MMCS scattering in a disc perpendicular to the shower axis
122 if (f[246]==0 && f[247]>0)
123 fImpactMax = f[247];
124
125 // CORSIKA scattering in a rectangle on the ground
126 if (f[246]>0 && f[247]>0)
127 *fLog << warn << "WARNING - Events scattered in a rectangle on the ground." << endl;
128
129 // Implemented in CORSIKA Version >= 6.822
130 memcpy(fAtmosphericLayers, f+248, 5*4);
131
132 memcpy(fAtmosphericCoeffA, f+253, 5*4);
133 memcpy(fAtmosphericCoeffB, f+258, 5*4);
134 memcpy(fAtmosphericCoeffC, f+263, 5*4);
135
136 // -------------------- Read first event header -------------------
137
138 // FIXME: Add sanity checks!
139
140 // f[76] Cherenkov flag:
141 // bit(1) : CERENKOV option compiled in
142 // bit(2) : IACT option compiled in
143 // bit(3) : CEFFIC option compiled in
144 // bit(4) : ATMEXT option compiled in
145 // bit(5) : ATMEXT option used with refraction enabled
146 // bit(6) : VOLUMEDET option compiled in
147 // bit(7) : CURVED option compiled in
148 // bit(9) : SLATN option compiled in
149 // 11-21 : table number for externam athmosphere (but<1024)
150 //
151 // f[78] Curved athmosphere? (0=flat, 1=curved)
152 // f[84] cherenkov bunch size
153 // f[93] flag for additinal muon information of particle output file
154 // f[145] Muon multiple scattering flag
155
156 if (!fInFormat->SeekNextBlock("EVTH", 1202))
157 return kFALSE;
158
159 Float_t g[273];
160 if (!fInFormat->ReadData(272, g))
161 return kFALSE;
162
163 fInFormat->UnreadLastData();
164 fInFormat->UnreadLastHeader();
165
166 const Int_t n = TMath::Nint(g[96]); // Number i of uses of each cherenkov event
167 if (n!=1)
168 {
169 *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl;
170 *fLog << warn << "WARNING - This error is replaced by a warning." << endl;
171 }
172
173 fParticleID = TMath::Nint(g[1]);
174
175 // MAGNETIC FIELD: x/z-component of earth magnetic field in muT
176 fMagneticFieldX = g[69]; // x-component ( BX)
177 fMagneticFieldZ = -g[70]; // z-component (-BZ)
178 fMagneticFieldAz = g[91]; // Azimuth angle of magnetic north expressed in telescope coordinates
179
180 // WITH rounding: unbelievable!
181 fCerenkovFlag = TMath::Nint(g[75]);
182
183 fZdMin = g[79]; // lower edge of theta in °
184 fZdMax = g[80]; // upper edge of theta in °
185 fAzMin = 180-g[81]; // lower edge of phi in °
186 fAzMax = 180-g[82]; // upper edge of phi in °
187 // FIXME: Correct for direction of magnetic field!
188
189 if (TMath::Nint(g[83])!=1)
190 *fLog << warn << "WARNING - Cherenkov bunch size not 1, but " << g[83] << endl;
191
192 // g[84] Number of cherenkov detectors in x
193 // g[85] Number of cherenkov detectors in y
194 // g[86] Grid spacing x
195 // g[87] Grid spacing y
196 // g[88] Length of detectors in x
197 // g[89] Length of detectors in y
198
199 fImpactMax = -1;
200/*
201 // This is a trick to use CERARY for storage of the
202 // maximum simulated impact
203 if (TMath::Nint(g[84])==1 && TMath::Nint(g[85])==1 &&
204 TMath::Nint(g[88])==1 && TMath::Nint(g[89])==1 &&
205 g[86]==g[87])
206 fImpactMax = g[86];
207 */
208 fWavelengthMin = g[94]; // Cherenkov bandwidth lower end in nm
209 fWavelengthMax = g[95]; // Cherenkov bandwidth upper end in nm
210
211 fViewConeInnerAngle = g[151]; // inner angle of view cone (°)
212 fViewConeOuterAngle = g[152]; // outer angle of view cone (°)
213
214 return kTRUE;
215}
216
217Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
218{
219
220 if (!fInFormat->SeekNextBlock("RUNE", 1210))
221 return kFALSE;
222
223 Float_t f[2];
224 if (!fInFormat->ReadData(2, f))
225 return kFALSE;
226
227 const UInt_t runnum = TMath::Nint(f[0]);
228 if (runnum!=fRunNumber)
229 {
230 *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
231 *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
232 return kFALSE;
233 }
234
235 fNumEvents = TMath::Nint(f[1]);
236
237 return kTRUE;
238}
239
240Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin)
241{
242 // Search subblockwise backward (Block: 5733*4 = 21*273*4)
243 for (int i=1; i<22; i++)
244 {
245 fin.seekg(-i*273*4, ios::end);
246
247 char runh[4];
248 fin.read(runh, 4);
249
250 if (!memcmp(runh, "RUNE", 4))
251 {
252 fin.seekg(-4, ios::cur);
253 return kTRUE;
254 }
255 }
256
257 return kFALSE;
258}
259
260// --------------------------------------------------------------------------
261//
262// print run header information on *fLog. The option 'header' supresses
263// the pixel index translation table.
264//
265void MCorsikaRunHeader::Print(Option_t *t) const
266{
267 *fLog << all << endl;
268 *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
269 *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl;
270 if (fNumEvents>0)
271 *fLog << "Num Events: " << fNumEvents << endl;
272 *fLog << "Obs Level: ";
273 for (Byte_t i=0; i<fNumObsLevel; i++)
274 *fLog << " " << fObsLevel[i]/100. << "m";
275 *fLog << endl;
276
277 *fLog << "MagneticField: X/Z=(" << fMagneticFieldX << "/";
278 *fLog << fMagneticFieldZ << ")" << UTF8::kMu << "T Az=" << fMagneticFieldAz*TMath::RadToDeg() << UTF8::kDeg << " (magnetic North w.r.t. North)" << endl;
279
280 *fLog << "Spectrum: Slope=" << fSlopeSpectrum << " (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" << endl;
281 *fLog << "Wavelength: " << fWavelengthMin << "nm - " << fWavelengthMax << "nm" << endl;
282
283 if (fImpactMax>0)
284 *fLog << "ImpactMax: " << fImpactMax << "cm" << endl;
285 if (fViewConeOuterAngle>0)
286 *fLog << "ViewCone: " << fViewConeInnerAngle << UTF8::kDeg << " - " << fViewConeOuterAngle << UTF8::kDeg << endl;
287
288 if (fZdMax>=0)
289 {
290 *fLog << "Zd/Az: " << fZdMin << UTF8::kDeg;
291 if (fZdMin==fZdMax)
292 *fLog << " (fixed)";
293 else
294 *fLog << "-" << fZdMax << UTF8::kDeg;
295 *fLog << " / " << fAzMin << UTF8::kDeg;
296 if (fAzMin==fAzMax)
297 *fLog << " (fixed)";
298 else
299 *fLog << "-" << fAzMax << UTF8::kDeg;
300 *fLog << " w.r.t. magnetic North." << endl;
301 }
302
303 if (fImpactMax>0)
304 *fLog << "Max.sim.Impact: " << fImpactMax << "cm" << endl;
305
306 *fLog << "Options used: ";
307 if (Has(kCerenkov))
308 *fLog << " CERENKOV";
309 if (Has(kIact))
310 *fLog << " IACT";
311 if (Has(kCeffic))
312 *fLog << " CEFFIC";
313 if (Has(kAtmext))
314 *fLog << " ATMEXT" << GetNumAtmosphericModel();
315 if (Has(kRefraction))
316 *fLog << " +Refraction";
317 if (Has(kVolumedet))
318 *fLog << " VOLUMEDET";
319 if (Has(kCurved))
320 *fLog << " CURVED";
321 if (Has(kSlant))
322 *fLog << " SLANT";
323 *fLog << endl;
324
325 if (HasLayers())
326 {
327 *fLog << "Atm.Layers: ";
328 for (int i=0; i<5; i++)
329 *fLog << " " << fAtmosphericLayers[i];
330 }
331 *fLog << "Atm.Coeff A: ";
332 for (int i=0; i<5; i++)
333 *fLog << " " << fAtmosphericCoeffA[i];
334 *fLog << endl;
335 *fLog << "Atm.Coeff B: ";
336 for (int i=0; i<5; i++)
337 *fLog << " " << fAtmosphericCoeffB[i];
338 *fLog << endl;
339 *fLog << "Atm.Coeff C: ";
340 for (int i=0; i<5; i++)
341 *fLog << " " << fAtmosphericCoeffC[i];
342 *fLog << endl;
343
344}
345
Note: See TracBrowser for help on using the repository browser.