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

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