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

Last change on this file since 20095 was 20044, checked in by tbretz, 4 years ago
Added telescope number, increased version accordingly.
File size: 13.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// Class Version 3:
46// ----------------
47// + UInt_t fNumReuse
48//
49// Class Version 4:
50// ----------------
51// + UInt_t fCerenkovFileOption
52// + UInt_t fHadronModelLowEnergy
53// + UInt_t fHadronModelHighEnergy
54// + Float_t fTransitionEnergy
55// + Bool_t fCurvedAtmosphere
56// + Float_t fEnergyCutoffHadrons
57// + Float_t fEnergyCutoffMuons
58// + Float_t fEnergyCutoffElectrons
59// + Float_t fEnergyCutoffPhotons
60// + Float_t fThinningEnergyFractionH
61// + Float_t fThinningEnergyFractionEM
62// + Float_t fThinningWeightLimitH
63// + Float_t fThinningWeightLimitEM
64// + Float_t fThinningMaxRadius
65//
66// Class Version 5:
67// ----------------
68// + fTelescopeNumber
69//
70////////////////////////////////////////////////////////////////////////////
71
72#include "MCorsikaRunHeader.h"
73#include "MCorsikaFormat.h"
74
75#include <fstream>
76#include <iomanip>
77
78#include "MLog.h"
79#include "MLogManip.h"
80
81#include "MMcEvt.hxx"
82
83ClassImp(MCorsikaRunHeader);
84
85using namespace std;
86
87const Double_t MCorsikaRunHeader::fgEarthRadius = 637131500; // [cm] Earth radius as defined in CORSIKA
88
89// --------------------------------------------------------------------------
90//
91// Default constructor. Creates array which stores the pixel assignment.
92//
93//
94MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title)
95 : fNumObsLevel(0), fImpactMax(-1), fZdMin(0), fZdMax(-1),
96 fAzMin(0), fAzMax(0), fWavelengthMin(-1), fWavelengthMax(-1),
97 fViewConeInnerAngle(0), fViewConeOuterAngle(-1)
98{
99 fName = name ? name : "MCorsikaRunHeader";
100 fTitle = title ? title : "Raw Run Header Information";
101}
102
103// --------------------------------------------------------------------------
104//
105// Read in one run header from the binary file
106//
107Bool_t MCorsikaRunHeader::ReadEvt(MCorsikaFormat * fInFormat, const uint32_t &blockLength)
108{
109 vector<Float_t> f(blockLength);
110 if (!fInFormat->Read(f.data(), blockLength*sizeof(Float_t)))
111 return kFALSE;
112
113 fTelescopeNumber = 0;
114 fRunNumber = TMath::Nint(f[0]);
115 fNumEvents = 0;
116
117 fRunStart.SetCorsikaTime(f[1]);
118
119 fProgramVersion = f[2];
120 fNumObsLevel = TMath::Nint(f[3]);
121
122 if (fNumObsLevel!=1)
123 {
124 *fLog << err << "ERROR - Currently only one observation level is allowed." << endl;
125 return kFALSE;
126 }
127
128 memset(fObsLevel, 0, 10*4);
129 memcpy(fObsLevel, f.data()+4, fNumObsLevel*4);
130
131 fSlopeSpectrum = f[14];
132 fEnergyMin = f[15];
133 fEnergyMax = f[16];
134
135 fEnergyCutoffHadrons = f[17];
136 fEnergyCutoffMuons = f[18];
137 fEnergyCutoffElectrons = f[19];
138 fEnergyCutoffPhotons = f[20];
139
140 // Implemented in CORSIKA Version >= 6.822
141 fImpactMax = -1;
142
143 // CORSIKA scattering in a disc on the ground
144 if (f[246]>0 && f[247]==0 && !fInFormat->IsEventioFormat())
145 {
146 *fLog << warn << "WARNING - Events scattered in a disc on the ground." << endl;
147 fImpactMax = f[246];
148 }
149
150 // MMCS scattering in a disc perpendicular to the shower axis
151 if (f[246]==0 && f[247]>0)
152 fImpactMax = f[247];
153
154 // CORSIKA scattering in a rectangle on the ground
155 if (f[246]>0 && f[247]>0)
156 *fLog << warn << "WARNING - Events scattered in a rectangle on the ground." << endl;
157
158 // Implemented in CORSIKA Version >= 6.822
159 memcpy(fAtmosphericLayers, f.data()+248, 5*4);
160
161 memcpy(fAtmosphericCoeffA, f.data()+253, 5*4);
162 memcpy(fAtmosphericCoeffB, f.data()+258, 5*4);
163 memcpy(fAtmosphericCoeffC, f.data()+263, 5*4);
164
165 return kTRUE;
166}
167
168// --------------------------------------------------------------------------
169//
170// Read in one event header. It is called for the first event header after
171// a run header
172//
173Bool_t MCorsikaRunHeader::ReadEventHeader(Float_t * g)
174{
175
176 // -------------------- Read first event header -------------------
177
178 // FIXME: Add sanity checks!
179
180 // f[76] Cherenkov flag:
181 // bit(1) : CERENKOV option compiled in
182 // bit(2) : IACT option compiled in
183 // bit(3) : CEFFIC option compiled in
184 // bit(4) : ATMEXT option compiled in
185 // bit(5) : ATMEXT option used with refraction enabled
186 // bit(6) : VOLUMEDET option compiled in
187 // bit(7) : CURVED option compiled in
188 // bit(9) : SLATN option compiled in
189 // 11-21 : table number for externam athmosphere (but<1024)
190 //
191 // f[78] Curved athmosphere? (0=flat, 1=curved)
192 // f[84] cherenkov bunch size
193 // f[93] flag for additinal muon information of particle output file
194 // f[145] Muon multiple scattering flag
195
196
197 fNumReuse = TMath::Nint(g[96]); // Number i of uses of each cherenkov event
198
199 fParticleID = TMath::Nint(g[1]);
200
201 // MAGNETIC FIELD: x/z-component of earth magnetic field in muT
202 fMagneticFieldX = g[69]; // x-component ( BX)
203 fMagneticFieldZ = -g[70]; // z-component (-BZ)
204 fMagneticFieldAz = g[91]; // Azimuth angle of magnetic north expressed in telescope coordinates
205
206 fHadronModelLowEnergy = TMath::Nint(g[73]);
207 fHadronModelHighEnergy = TMath::Nint(g[74]);
208
209 // WITH rounding: unbelievable!
210 fCerenkovFlag = TMath::Nint(g[75]);
211 fCerenkovFileOption = TMath::Nint(g[90]);
212
213 fCurvedAtmosphere = TMath::Nint(g[77])==2;
214
215 fZdMin = g[79]; // lower edge of theta in °
216 fZdMax = g[80]; // upper edge of theta in °
217 fAzMin = 180-g[81]; // lower edge of phi in °
218 fAzMax = 180-g[82]; // upper edge of phi in °
219 // FIXME: Correct for direction of magnetic field!
220
221 if (TMath::Nint(g[83])!=1)
222 *fLog << warn << "WARNING - Cherenkov bunch size not 1, but " << g[83] << endl;
223
224 // g[84] Number of cherenkov detectors in x
225 // g[85] Number of cherenkov detectors in y
226 // g[86] Grid spacing x
227 // g[87] Grid spacing y
228 // g[88] Length of detectors in x
229 // g[89] Length of detectors in y
230
231 fImpactMax = -1;
232/*
233 // This is a trick to use CERARY for storage of the
234 // maximum simulated impact
235 if (TMath::Nint(g[84])==1 && TMath::Nint(g[85])==1 &&
236 TMath::Nint(g[88])==1 && TMath::Nint(g[89])==1 &&
237 g[86]==g[87])
238 fImpactMax = g[86];
239 */
240 fWavelengthMin = g[94]; // Cherenkov bandwidth lower end in nm
241 fWavelengthMax = g[95]; // Cherenkov bandwidth upper end in nm
242
243 fThinningEnergyFractionH = g[146]; // EFRCTHN
244 fThinningEnergyFractionEM = g[147]; // EFRCTHN*THINRAT
245 fThinningWeightLimitH = g[148]; // WMAX
246 fThinningWeightLimitEM = g[149]; // WMAX*WEITRAT
247 fThinningMaxRadius = g[150]; // Max radial radius for thinning
248
249 fViewConeInnerAngle = g[151]; // inner angle of view cone (°)
250 fViewConeOuterAngle = g[152]; // outer angle of view cone (°)
251
252 fTransitionEnergy = g[153];
253
254 return kTRUE;
255}
256
257Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat, Bool_t runNumberVerify)
258{
259 Float_t f[2];
260 if (!fInFormat->Read(f, 2 * sizeof(Float_t)))
261 return kFALSE;
262
263 if (runNumberVerify)
264 {
265 const UInt_t runnum = TMath::Nint(f[0]);
266 if (runnum!=fRunNumber)
267 {
268 *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
269 *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
270 return kFALSE;
271 }
272 }
273
274 fNumEvents = TMath::Nint(f[1]);
275
276 return kTRUE;
277}
278
279// --------------------------------------------------------------------------
280//
281// print run header information on *fLog. The option 'header' supresses
282// the pixel index translation table.
283//
284void MCorsikaRunHeader::Print(Option_t *t) const
285{
286 *fLog << all << endl;
287 *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
288 *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl;
289 if (fNumEvents>0)
290 *fLog << "Num Events: " << fNumEvents << " (reuse " << fNumReuse << " times)" << endl;
291 *fLog << "Obs Level: ";
292 for (Byte_t i=0; i<fNumObsLevel; i++)
293 *fLog << " " << fObsLevel[i]/100. << "m";
294 *fLog << endl;
295
296 *fLog << "MagneticField: X/Z=(" << fMagneticFieldX << "/";
297 *fLog << fMagneticFieldZ << ")" << UTF8::kMu << "T Az=" << fMagneticFieldAz*TMath::RadToDeg() << UTF8::kDeg << " (magnetic North w.r.t. North)" << endl;
298
299 *fLog << "Spectrum: Slope=" << fSlopeSpectrum << " (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" << endl;
300 *fLog << "Wavelength: " << fWavelengthMin << "nm - " << fWavelengthMax << "nm" << endl;
301
302 if (fImpactMax>0)
303 *fLog << "ImpactMax: " << fImpactMax << "cm" << endl;
304 if (fViewConeOuterAngle>0)
305 *fLog << "ViewCone: " << fViewConeInnerAngle << UTF8::kDeg << " - " << fViewConeOuterAngle << UTF8::kDeg << endl;
306
307 *fLog << "Zd/Az: ";
308 if (fZdMax>=0 && fZdMin<360)
309 {
310 *fLog << fZdMin << UTF8::kDeg;
311 if (fZdMin==fZdMax)
312 *fLog << " (fixed)";
313 else
314 *fLog << "-" << fZdMax << UTF8::kDeg;
315 *fLog << " / " << fAzMin << UTF8::kDeg;
316 if (fAzMin==fAzMax)
317 *fLog << " (fixed)";
318 else
319 *fLog << "-" << fAzMax << UTF8::kDeg;
320 *fLog << " w.r.t. magnetic North." << endl;
321 }
322
323 if (fZdMin>=360) // 4010.7
324 *fLog << "-trajectory-" << endl;
325
326
327 if (fImpactMax>0)
328 *fLog << "Max.sim.Impact: " << fImpactMax << "cm" << endl;
329
330 *fLog << "Energy cutoff: ";
331 *fLog << fEnergyCutoffHadrons << "GeV (hadrons), ";
332 *fLog << fEnergyCutoffMuons << "GeV (muons), ";
333 *fLog << fEnergyCutoffElectrons << "GeV (electrons), ";
334 *fLog << fEnergyCutoffPhotons << "GeV (photons)";
335 *fLog << endl;
336
337 *fLog << "Thinning: ";
338 if (fThinningWeightLimitH>0)
339 {
340 *fLog << "HADRONIC: E/Eth>" << fThinningEnergyFractionH << " (w>" << fThinningWeightLimitH << "), ";
341 *fLog << "EM: E/Eth>" << fThinningEnergyFractionEM << " (w>" << fThinningWeightLimitEM << "), ";
342 *fLog << "R>" << fThinningMaxRadius << "cm";
343 *fLog << endl;
344 }
345 else
346 *fLog << "<off>" << endl;
347
348 *fLog << "Interact.model: ";
349 switch (fHadronModelLowEnergy)
350 {
351 case 1: *fLog << "GEISHA"; break;
352 case 2: *fLog << "UrQMD"; break;
353 case 3: *fLog << "FLUKA"; break;
354 default: *fLog << "<n/a>"; break;
355 }
356 *fLog << " / ";
357 switch (fHadronModelHighEnergy)
358 {
359 case 0: *fLog << "HDPM"; break;
360 case 1: *fLog << "VENUS"; break;
361 case 2: *fLog << "SIBYLL"; break;
362 case 3: *fLog << "QGSJET"; break;
363 case 4: *fLog << "DPMJET"; break;
364 case 5: *fLog << "neXus"; break;
365 case 6: *fLog << "EPOS"; break;
366 default: *fLog << "<n/a>"; break;
367 }
368 *fLog << " [lo/hi], Transition at " << fTransitionEnergy << " GeV" << endl;
369
370 *fLog << "Options used: ";
371 if (Has(kCerenkov))
372 *fLog << " CERENKOV";
373 if (Has(kIact))
374 *fLog << " IACT";
375 if (Has(kCeffic))
376 *fLog << " CEFFIC";
377 if (Has(kAtmext))
378 *fLog << " ATMEXT" << GetNumAtmosphericModel();
379 if (Has(kRefraction))
380 *fLog << " +Refraction";
381 if (Has(kVolumedet))
382 *fLog << " VOLUMEDET";
383 if (Has(kCurved))
384 *fLog << " CURVED" << (fCurvedAtmosphere?"<on>":"<off>");
385 if (Has(kSlant))
386 *fLog << " SLANT";
387 *fLog << " [" << hex << fCerenkovFlag << "]" << dec << endl;
388
389 if (Has(kCerenkov))
390 {
391 *fLog << "File format: ";
392 switch (fCerenkovFileOption)
393 {
394 case 0:
395 *fLog << "Cerenkov photons written to DAT-file.";
396 break;
397 case 1:
398 *fLog << "Cerenkov photons written to CER-file";
399 break;
400 case 2:
401 *fLog << "Cerenkov photons written to CER-file / Wavelength as 8th item in THIN option";
402 break;
403 default:
404 *fLog << "Cerenkov photons written to CER-file / Prod. height replaced by distance to array center.";
405 break;
406 }
407 *fLog << " [MCERFI=" << fCerenkovFileOption << "]" << endl;
408 }
409
410 if (HasLayers())
411 {
412 *fLog << "Atm.Layers: ";
413 for (int i=0; i<5; i++)
414 *fLog << " " << fAtmosphericLayers[i];
415 }
416 *fLog << endl;
417 *fLog << "Atm.Coeff A: ";
418 for (int i=0; i<5; i++)
419 *fLog << " " << fAtmosphericCoeffA[i];
420 *fLog << endl;
421 *fLog << "Atm.Coeff B: ";
422 for (int i=0; i<5; i++)
423 *fLog << " " << fAtmosphericCoeffB[i];
424 *fLog << endl;
425 *fLog << "Atm.Coeff C: ";
426 for (int i=0; i<5; i++)
427 *fLog << " " << fAtmosphericCoeffC[i];
428 *fLog << endl;
429
430
431}
432
Note: See TracBrowser for help on using the repository browser.