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

Last change on this file since 9942 was 9937, checked in by tbretz, 15 years ago
Added support for re-use of showers in MMCS.
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))
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))
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
224 if (!fInFormat->SeekNextBlock("RUNE", 1210))
225 return kFALSE;
226
227 Float_t f[2];
228 if (!fInFormat->ReadData(2, f))
229 return kFALSE;
230
231 const UInt_t runnum = TMath::Nint(f[0]);
232 if (runnum!=fRunNumber)
233 {
234 *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
235 *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
236 return kFALSE;
237 }
238
239 fNumEvents = TMath::Nint(f[1]);
240
241 return kTRUE;
242}
243
244Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin)
245{
246 // Search subblockwise backward (Block: 5733*4 = 21*273*4)
247 for (int i=1; i<22; i++)
248 {
249 fin.seekg(-i*273*4, ios::end);
250
251 char runh[4];
252 fin.read(runh, 4);
253
254 if (!memcmp(runh, "RUNE", 4))
255 {
256 fin.seekg(-4, ios::cur);
257 return kTRUE;
258 }
259 }
260
261 return kFALSE;
262}
263
264// --------------------------------------------------------------------------
265//
266// print run header information on *fLog. The option 'header' supresses
267// the pixel index translation table.
268//
269void MCorsikaRunHeader::Print(Option_t *t) const
270{
271 *fLog << all << endl;
272 *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
273 *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl;
274 if (fNumEvents>0)
275 *fLog << "Num Events: " << fNumEvents << " (reuse " << fNumReuse << " times)" << endl;
276 *fLog << "Obs Level: ";
277 for (Byte_t i=0; i<fNumObsLevel; i++)
278 *fLog << " " << fObsLevel[i]/100. << "m";
279 *fLog << endl;
280
281 *fLog << "MagneticField: X/Z=(" << fMagneticFieldX << "/";
282 *fLog << fMagneticFieldZ << ")" << UTF8::kMu << "T Az=" << fMagneticFieldAz*TMath::RadToDeg() << UTF8::kDeg << " (magnetic North w.r.t. North)" << endl;
283
284 *fLog << "Spectrum: Slope=" << fSlopeSpectrum << " (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" << endl;
285 *fLog << "Wavelength: " << fWavelengthMin << "nm - " << fWavelengthMax << "nm" << endl;
286
287 if (fImpactMax>0)
288 *fLog << "ImpactMax: " << fImpactMax << "cm" << endl;
289 if (fViewConeOuterAngle>0)
290 *fLog << "ViewCone: " << fViewConeInnerAngle << UTF8::kDeg << " - " << fViewConeOuterAngle << UTF8::kDeg << endl;
291
292 if (fZdMax>=0)
293 {
294 *fLog << "Zd/Az: " << fZdMin << UTF8::kDeg;
295 if (fZdMin==fZdMax)
296 *fLog << " (fixed)";
297 else
298 *fLog << "-" << fZdMax << UTF8::kDeg;
299 *fLog << " / " << fAzMin << UTF8::kDeg;
300 if (fAzMin==fAzMax)
301 *fLog << " (fixed)";
302 else
303 *fLog << "-" << fAzMax << UTF8::kDeg;
304 *fLog << " w.r.t. magnetic North." << endl;
305 }
306
307 if (fImpactMax>0)
308 *fLog << "Max.sim.Impact: " << fImpactMax << "cm" << endl;
309
310 *fLog << "Options used: ";
311 if (Has(kCerenkov))
312 *fLog << " CERENKOV";
313 if (Has(kIact))
314 *fLog << " IACT";
315 if (Has(kCeffic))
316 *fLog << " CEFFIC";
317 if (Has(kAtmext))
318 *fLog << " ATMEXT" << GetNumAtmosphericModel();
319 if (Has(kRefraction))
320 *fLog << " +Refraction";
321 if (Has(kVolumedet))
322 *fLog << " VOLUMEDET";
323 if (Has(kCurved))
324 *fLog << " CURVED";
325 if (Has(kSlant))
326 *fLog << " SLANT";
327 *fLog << endl;
328
329 if (HasLayers())
330 {
331 *fLog << "Atm.Layers: ";
332 for (int i=0; i<5; i++)
333 *fLog << " " << fAtmosphericLayers[i];
334 }
335 *fLog << "Atm.Coeff A: ";
336 for (int i=0; i<5; i++)
337 *fLog << " " << fAtmosphericCoeffA[i];
338 *fLog << endl;
339 *fLog << "Atm.Coeff B: ";
340 for (int i=0; i<5; i++)
341 *fLog << " " << fAtmosphericCoeffB[i];
342 *fLog << endl;
343 *fLog << "Atm.Coeff C: ";
344 for (int i=0; i<5; i++)
345 *fLog << " " << fAtmosphericCoeffC[i];
346 *fLog << endl;
347
348}
349
Note: See TracBrowser for help on using the repository browser.