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

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