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

Last change on this file since 9352 was 9352, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 9.5 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 fAtmosphericCoeffA[5]
40// + Float_t fAtmosphericCoeffB[5]
41// + Float_t fAtmosphericCoeffC[5]
42// + UInt_t fCerenkovFlag
43//
44////////////////////////////////////////////////////////////////////////////
45
46#include "MCorsikaRunHeader.h"
47
48#include <fstream>
49#include <iomanip>
50
51#include "MLog.h"
52#include "MLogManip.h"
53
54#include "MMcEvt.hxx"
55
56ClassImp(MCorsikaRunHeader);
57
58using namespace std;
59
60const Double_t MCorsikaRunHeader::fgEarthRadius = 637131500; // [cm] Earth radius as defined in CORSIKA
61
62// --------------------------------------------------------------------------
63//
64// Default constructor. Creates array which stores the pixel assignment.
65//
66//
67MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title)
68 : fNumObsLevel(0), fZdMin(0), fZdMax(-1), fAzMin(0), fAzMax(0),
69 fViewConeInnerAngle(0), fViewConeOuterAngle(-1)
70{
71 fName = name ? name : "MCorsikaRunHeader";
72 fTitle = title ? title : "Raw Run Header Information";
73}
74
75// --------------------------------------------------------------------------
76//
77// Read in one run header from the binary file
78//
79Bool_t MCorsikaRunHeader::ReadEvt(istream& fin)
80{
81 char runh[4];
82 fin.read(runh, 4);
83 if (memcmp(runh, "RUNH", 4))
84 {
85 *fLog << err << "ERROR - Wrong identifier: RUNH expected." << endl;
86 return kFALSE;
87 }
88
89 Float_t f[272*4];
90 fin.read((char*)f, 272*4);
91
92 fRunNumber = TMath::Nint(f[0]);
93 fNumEvents = 0;
94
95 fRunStart.SetCorsikaTime(f[1]);
96
97 fProgramVersion = f[2];
98 fNumObsLevel = TMath::Nint(f[3]);
99
100 if (fNumObsLevel!=1)
101 {
102 *fLog << err << "ERROR - Currently only one observation level is allowed." << endl;
103 return kFALSE;
104 }
105
106 memset(fObsLevel, 0, 10*4);
107 memcpy(fObsLevel, f+4, fNumObsLevel*4);
108
109 fSlopeSpectrum = f[14];
110 fEnergyMin = f[15];
111 fEnergyMax = f[16];
112
113 // 0/10 not supported?
114
115 memcpy(fAtmosphericCoeffA, f+253, 5*4);
116 memcpy(fAtmosphericCoeffB, f+258, 5*4);
117 memcpy(fAtmosphericCoeffC, f+263, 5*4);
118
119 // -------------------- Read first event header -------------------
120
121 // FIXME: Add sanity checks!
122
123 // f[76] Cherenkov flag:
124 // bit(1) : CERENKOV option compiled in
125 // bit(2) : IACT option compiled in
126 // bit(3) : CEFFIC option compiled in
127 // bit(4) : ATMEXT option compiled in
128 // bit(5) : ATMEXT option used with refraction enabled
129 // bit(6) : VOLUMEDET option compiled in
130 // bit(7) : CURVED option compiled in
131 // bit(9) : SLATN option compiled in
132 // 11-21 : table number for externam athmosphere (but<1024)
133 //
134 // f[78] Curved athmosphere? (0=flat, 1=curved)
135 // f[84] cherenkov bunch size
136 // f[93] flag for additinal muon information of particle output file
137 // f[145] Muon multiple scattering flag
138 // f[156] altitude of horizontal shower axis
139
140 char evth[4];
141 fin.read(evth, 4);
142 if (memcmp(evth, "EVTH", 4))
143 {
144 *fLog << err << "ERROR - Wrong identifier: EVTH expected." << endl;
145 return kFALSE;
146 }
147
148 Float_t g[273];
149 fin.read((char*)&g, 273*4);
150 if (fin.eof())
151 return kFALSE;
152
153 fin.seekg(-274*4, ios::cur);
154
155 const Int_t n = TMath::Nint(g[96]); // Number i of uses of each cherenkov event
156 if (n!=1)
157 {
158 *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl;
159 return kFALSE;
160 }
161
162 fParticleID = TMath::Nint(g[1]);
163
164 // MAGNETIC FIELD: x/z-component of earth magnetic field in muT
165 fMagneticFieldX = g[69]; // x-component ( BX)
166 fMagneticFieldZ = -g[70]; // z-component (-BZ)
167 fMagneticFieldAz = g[91]; // Azimuth angle of magnetic north expressed in telescope coordinates
168
169 // WITH rounding: unbelievable!
170 fCerenkovFlag = TMath::Nint(g[75]);
171
172 fZdMin = g[79]; // lower edge of theta in °
173 fZdMax = g[80]; // upper edge of theta in °
174 fAzMin = 180-g[81]; // lower edge of phi in °
175 fAzMax = 180-g[82]; // upper edge of phi in °
176 // FIXME: Correct for direction of magnetic field!
177
178 // Number of cherenkov detectors in x
179 // Number of cherenkov detectors in y
180 // Grid spacing x
181 // Grid spacing y
182
183 // g[93] angle between array x-direction and magnetic north
184
185
186 fImpactMax = g[86];
187
188 fWavelengthMin = g[94]; // Cherenkov bandwidth lower end in nm
189 fWavelengthMax = g[95]; // Cherenkov bandwidth upper end in nm
190
191 fViewConeInnerAngle = g[151]; // inner angle of view cone (°)
192 fViewConeOuterAngle = g[152]; // outer angle of view cone (°)
193
194 return kTRUE;
195}
196
197Bool_t MCorsikaRunHeader::ReadEvtEnd(istream& fin)
198{
199 char runh[4];
200 fin.read(runh, 4);
201 if (memcmp(runh, "RUNE", 4))
202 {
203 *fLog << err << "ERROR - Wrong identifier: RUNE expected." << endl;
204 return kFALSE;
205 }
206
207 Float_t f[2];
208 fin.read((char*)f, 2*4);
209
210 const UInt_t runnum = TMath::Nint(f[0]);
211 if (runnum!=fRunNumber)
212 {
213 *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
214 *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
215 return kFALSE;
216 }
217
218 fNumEvents = TMath::Nint(f[1]);
219
220 fin.seekg(270*4, ios::cur); // skip the remaining data of this block
221
222 return kTRUE;
223}
224
225Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin)
226{
227 // Search subblockwise backward (Block: 5733*4 = 21*273*4)
228 for (int i=1; i<22; i++)
229 {
230 fin.seekg(-i*273*4, ios::end);
231
232 char runh[4];
233 fin.read(runh, 4);
234
235 if (!memcmp(runh, "RUNE", 4))
236 {
237 fin.seekg(-4, ios::cur);
238 return kTRUE;
239 }
240 }
241
242 return kFALSE;
243}
244
245// --------------------------------------------------------------------------
246//
247// print run header information on *fLog. The option 'header' supresses
248// the pixel index translation table.
249//
250void MCorsikaRunHeader::Print(Option_t *t) const
251{
252 *fLog << all << endl;
253 *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
254 *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl;
255 if (fNumEvents>0)
256 *fLog << "Num Events: " << fNumEvents << endl;
257 *fLog << "Obs Level: ";
258 for (Byte_t i=0; i<fNumObsLevel; i++)
259 *fLog << " " << fObsLevel[i]/100. << "m";
260 *fLog << endl;
261
262 *fLog << "MagneticField: X/Z=(" << fMagneticFieldX << "/";
263 *fLog << fMagneticFieldZ << ")" << UTF8::kMu << "T Az=" << fMagneticFieldAz*TMath::RadToDeg() << "deg (magnetic North w.r.t. North)" << endl;
264
265 *fLog << "Spectrum: Slope=" << fSlopeSpectrum << " (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" << endl;
266 *fLog << "Wavelength: " << fWavelengthMin << "nm - " << fWavelengthMax << "nm" << endl;
267
268 if (fViewConeOuterAngle>0)
269 *fLog << "ViewCone: " << fViewConeInnerAngle << UTF8::kDeg << "-" << fViewConeOuterAngle << UTF8::kDeg << endl;
270
271 if (fZdMax>=0)
272 {
273 *fLog << "Zd/Az: " << fZdMin << UTF8::kDeg;
274 if (fZdMin==fZdMax)
275 *fLog << " (fixed)";
276 else
277 *fLog << "-" << fZdMax << UTF8::kDeg;
278 *fLog << " / " << fAzMin << UTF8::kDeg;
279 if (fAzMin==fAzMax)
280 *fLog << " (fixed)";
281 else
282 *fLog << "-" << fAzMax << UTF8::kDeg;
283 *fLog << " w.r.t. magnetic North." << endl;
284 }
285
286 *fLog << "Options used: ";
287 if (Has(kCerenkov))
288 *fLog << " CERENKOV";
289 if (Has(kIact))
290 *fLog << " IACT";
291 if (Has(kCeffic))
292 *fLog << " CEFFIC";
293 if (Has(kAtmext))
294 *fLog << " ATMEXT";
295 if (Has(kRefraction))
296 *fLog << " AtmRefraction";
297 if (Has(kVolumedet))
298 *fLog << " VOLUMEDET";
299 if (Has(kCurved))
300 *fLog << " CURVED";
301 if (Has(kSlant))
302 *fLog << " SLANT";
303 *fLog << endl;
304
305 *fLog << "Atm.Coeff A: ";
306 for (int i=0; i<5; i++)
307 *fLog << " " << fAtmosphericCoeffA[i];
308 *fLog << endl;
309 *fLog << "Atm.Coeff B: ";
310 for (int i=0; i<5; i++)
311 *fLog << " " << fAtmosphericCoeffB[i];
312 *fLog << endl;
313 *fLog << "Atm.Coeff C: ";
314 for (int i=0; i<5; i++)
315 *fLog << " " << fAtmosphericCoeffC[i];
316 *fLog << endl;
317
318}
319
Note: See TracBrowser for help on using the repository browser.