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