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

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