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

Last change on this file since 9362 was 9362, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 9.9 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), fImpactMax(-1), fZdMin(0), fZdMax(-1),
69 fAzMin(0), fAzMax(0), 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 fImpactMax = -1;
186/*
187 // This is a trick to use CERARY for storage of the
188 // maximum simulated impact
189 if (TMath::Nint(g[84])==1 && TMath::Nint(g[85])==1 &&
190 TMath::Nint(g[88])==1 && TMath::Nint(g[89])==1 &&
191 g[86]==g[87])
192 fImpactMax = g[86];
193 */
194 fWavelengthMin = g[94]; // Cherenkov bandwidth lower end in nm
195 fWavelengthMax = g[95]; // Cherenkov bandwidth upper end in nm
196
197 fViewConeInnerAngle = g[151]; // inner angle of view cone (°)
198 fViewConeOuterAngle = g[152]; // outer angle of view cone (°)
199
200 return kTRUE;
201}
202
203Bool_t MCorsikaRunHeader::ReadEvtEnd(istream& fin)
204{
205 char runh[4];
206 fin.read(runh, 4);
207 if (memcmp(runh, "RUNE", 4))
208 {
209 *fLog << err << "ERROR - Wrong identifier: RUNE expected." << endl;
210 return kFALSE;
211 }
212
213 Float_t f[2];
214 fin.read((char*)f, 2*4);
215
216 const UInt_t runnum = TMath::Nint(f[0]);
217 if (runnum!=fRunNumber)
218 {
219 *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
220 *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
221 return kFALSE;
222 }
223
224 fNumEvents = TMath::Nint(f[1]);
225
226 fin.seekg(270*4, ios::cur); // skip the remaining data of this block
227
228 return kTRUE;
229}
230
231Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin)
232{
233 // Search subblockwise backward (Block: 5733*4 = 21*273*4)
234 for (int i=1; i<22; i++)
235 {
236 fin.seekg(-i*273*4, ios::end);
237
238 char runh[4];
239 fin.read(runh, 4);
240
241 if (!memcmp(runh, "RUNE", 4))
242 {
243 fin.seekg(-4, ios::cur);
244 return kTRUE;
245 }
246 }
247
248 return kFALSE;
249}
250
251// --------------------------------------------------------------------------
252//
253// print run header information on *fLog. The option 'header' supresses
254// the pixel index translation table.
255//
256void MCorsikaRunHeader::Print(Option_t *t) const
257{
258 *fLog << all << endl;
259 *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
260 *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl;
261 if (fNumEvents>0)
262 *fLog << "Num Events: " << fNumEvents << endl;
263 *fLog << "Obs Level: ";
264 for (Byte_t i=0; i<fNumObsLevel; i++)
265 *fLog << " " << fObsLevel[i]/100. << "m";
266 *fLog << endl;
267
268 *fLog << "MagneticField: X/Z=(" << fMagneticFieldX << "/";
269 *fLog << fMagneticFieldZ << ")" << UTF8::kMu << "T Az=" << fMagneticFieldAz*TMath::RadToDeg() << "deg (magnetic North w.r.t. North)" << endl;
270
271 *fLog << "Spectrum: Slope=" << fSlopeSpectrum << " (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" << endl;
272 *fLog << "Wavelength: " << fWavelengthMin << "nm - " << fWavelengthMax << "nm" << endl;
273
274 if (fViewConeOuterAngle>0)
275 *fLog << "ViewCone: " << fViewConeInnerAngle << UTF8::kDeg << " - " << fViewConeOuterAngle << UTF8::kDeg << endl;
276
277 if (fZdMax>=0)
278 {
279 *fLog << "Zd/Az: " << fZdMin << UTF8::kDeg;
280 if (fZdMin==fZdMax)
281 *fLog << " (fixed)";
282 else
283 *fLog << "-" << fZdMax << UTF8::kDeg;
284 *fLog << " / " << fAzMin << UTF8::kDeg;
285 if (fAzMin==fAzMax)
286 *fLog << " (fixed)";
287 else
288 *fLog << "-" << fAzMax << UTF8::kDeg;
289 *fLog << " w.r.t. magnetic North." << endl;
290 }
291
292 if (fImpactMax>0)
293 *fLog << "Max.sim.Impact: " << fImpactMax << "cm" << endl;
294
295 *fLog << "Options used: ";
296 if (Has(kCerenkov))
297 *fLog << " CERENKOV";
298 if (Has(kIact))
299 *fLog << " IACT";
300 if (Has(kCeffic))
301 *fLog << " CEFFIC";
302 if (Has(kAtmext))
303 *fLog << " ATMEXT" << GetNumAtmosphericModel();
304 if (Has(kRefraction))
305 *fLog << " AtmRefraction";
306 if (Has(kVolumedet))
307 *fLog << " VOLUMEDET";
308 if (Has(kCurved))
309 *fLog << " CURVED";
310 if (Has(kSlant))
311 *fLog << " SLANT";
312 *fLog << endl;
313
314 *fLog << "Atm.Coeff A: ";
315 for (int i=0; i<5; i++)
316 *fLog << " " << fAtmosphericCoeffA[i];
317 *fLog << endl;
318 *fLog << "Atm.Coeff B: ";
319 for (int i=0; i<5; i++)
320 *fLog << " " << fAtmosphericCoeffB[i];
321 *fLog << endl;
322 *fLog << "Atm.Coeff C: ";
323 for (int i=0; i<5; i++)
324 *fLog << " " << fAtmosphericCoeffC[i];
325 *fLog << endl;
326
327}
328
Note: See TracBrowser for help on using the repository browser.