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

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