source: trunk/Mars/mcorsika/MCorsikaRunHeader.cc@ 10009

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