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

Last change on this file since 19341 was 19341, checked in by tbretz, 8 months ago
Block length in number fo floats.
File size: 12.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// 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*sizeof(Float_t)))
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.