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

Last change on this file since 18546 was 18546, checked in by tbretz, 3 years ago
In EventIO the impact is not calculated on the ground, the trajectory option is detected differently, the cherenkov flag is output also as a hex number
File size: 10.5 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    Float_t f[272];
88    if (!fInFormat->Read(f, 272 * sizeof(Float_t)))
89        return kFALSE;
90
91    fRunNumber = TMath::Nint(f[0]);
92    fNumEvents = 0;
93
94    fRunStart.SetCorsikaTime(f[1]);
95
96    fProgramVersion = f[2];
97    fNumObsLevel    = TMath::Nint(f[3]);
98
99    if (fNumObsLevel!=1)
100    {
101        *fLog << err << "ERROR - Currently only one observation level is allowed." << endl;
102        return kFALSE;
103    }
104
105    memset(fObsLevel, 0, 10*4);
106    memcpy(fObsLevel, f+4, fNumObsLevel*4);
107
108    fSlopeSpectrum  = f[14];
109    fEnergyMin      = f[15];
110    fEnergyMax      = f[16];
111
112    // Implemented in CORSIKA Version >= 6.822
113    fImpactMax = -1;
114
115    // CORSIKA scattering in a disc on the ground
116    if (f[246]>0 && f[247]==0 && !fInFormat->IsEventioFormat())
117    {
118        *fLog << warn << "WARNING - Events scattered in a disc on the ground." << endl;
119        fImpactMax = f[246];
120    }
121
122    // MMCS scattering in a disc perpendicular to the shower axis
123    if (f[246]==0 && f[247]>0)
124        fImpactMax = f[247];
125
126    // CORSIKA scattering in a rectangle on the ground
127    if (f[246]>0 && f[247]>0)
128        *fLog << warn << "WARNING - Events scattered in a rectangle on the ground." << endl;
129
130    // Implemented in CORSIKA Version >= 6.822
131    memcpy(fAtmosphericLayers, f+248, 5*4);
132
133    memcpy(fAtmosphericCoeffA, f+253, 5*4);
134    memcpy(fAtmosphericCoeffB, f+258, 5*4);
135    memcpy(fAtmosphericCoeffC, f+263, 5*4);
136
137    return kTRUE;
138}
139
140// --------------------------------------------------------------------------
141//
142// Read in one event header. It is called for the first event header after 
143// a run header                                                             
144//
145Bool_t MCorsikaRunHeader::ReadEventHeader(Float_t * g)
146{
147
148    // -------------------- Read first event header -------------------
149
150    // FIXME: Add sanity checks!
151
152    // f[76] Cherenkov flag:
153    //        bit(1) : CERENKOV option compiled in
154    //        bit(2) : IACT option compiled in
155    //        bit(3) : CEFFIC option compiled in
156    //        bit(4) : ATMEXT option compiled in
157    //        bit(5) : ATMEXT option used with refraction enabled
158    //        bit(6) : VOLUMEDET option compiled in
159    //        bit(7) : CURVED option compiled in
160    //        bit(9) : SLATN option compiled in
161    //        11-21  : table number for externam athmosphere (but<1024)
162    //
163    // f[78]  Curved athmosphere? (0=flat, 1=curved)
164    // f[84]  cherenkov bunch size
165    // f[93]  flag for additinal muon information of particle output file
166    // f[145] Muon multiple scattering flag
167
168
169    fNumReuse = TMath::Nint(g[96]);  // Number i of uses of each cherenkov event
170
171    fParticleID = TMath::Nint(g[1]);
172
173    // MAGNETIC FIELD: x/z-component of earth magnetic field in muT
174    fMagneticFieldX  =  g[69];  // x-component ( BX)
175    fMagneticFieldZ  = -g[70];  // z-component (-BZ)
176    fMagneticFieldAz =  g[91];  // Azimuth angle of magnetic north expressed in telescope coordinates
177
178    // WITH rounding: unbelievable!
179    fCerenkovFlag = TMath::Nint(g[75]);
180
181    fZdMin = g[79];                // lower edge of theta in °
182    fZdMax = g[80];                // upper edge of theta in °
183    fAzMin = 180-g[81];            // lower edge of phi   in °
184    fAzMax = 180-g[82];            // upper edge of phi   in °
185    // FIXME: Correct for direction of magnetic field!
186
187    if (TMath::Nint(g[83])!=1)
188        *fLog << warn << "WARNING - Cherenkov bunch size not 1, but " << g[83] << endl;
189
190    // g[84] Number of cherenkov detectors in x
191    // g[85] Number of cherenkov detectors in y
192    // g[86] Grid spacing x
193    // g[87] Grid spacing y
194    // g[88] Length of detectors in x
195    // g[89] Length of detectors in y
196
197    fImpactMax = -1;
198/*
199    // This is a trick to use CERARY for storage of the
200    // maximum simulated impact
201    if (TMath::Nint(g[84])==1 && TMath::Nint(g[85])==1 &&
202        TMath::Nint(g[88])==1 && TMath::Nint(g[89])==1 &&
203        g[86]==g[87])
204        fImpactMax = g[86];
205 */
206    fWavelengthMin = g[94];        // Cherenkov bandwidth lower end in nm
207    fWavelengthMax = g[95];        // Cherenkov bandwidth upper end in nm
208
209    fViewConeInnerAngle = g[151];  // inner angle of view cone (°)
210    fViewConeOuterAngle = g[152];  // outer angle of view cone (°)
211
212    return kTRUE;
213}
214
215Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat, Bool_t runNumberVerify)
216{
217    Float_t f[2];
218    if (!fInFormat->Read(f, 2 * sizeof(Float_t)))
219        return kFALSE;
220
221    if (runNumberVerify)
222      {
223       const UInt_t runnum = TMath::Nint(f[0]);
224       if (runnum!=fRunNumber)
225       {
226           *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
227           *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
228           return kFALSE;
229       }
230      }
231
232    fNumEvents = TMath::Nint(f[1]);
233
234    return kTRUE;
235}
236
237// --------------------------------------------------------------------------
238//
239// print run header information on *fLog. The option 'header' supresses
240// the pixel index translation table.
241//
242void MCorsikaRunHeader::Print(Option_t *t) const
243{
244    *fLog << all << endl;
245    *fLog << "Run Number:     " << fRunNumber << "  (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
246    *fLog << "Particle ID:    " << MMcEvt::GetParticleName(fParticleID) << endl;
247    if (fNumEvents>0)
248        *fLog << "Num Events:     " << fNumEvents << " (reuse " << fNumReuse << " times)" << endl;
249    *fLog << "Obs Level:     ";
250    for (Byte_t i=0; i<fNumObsLevel; i++)
251        *fLog << " " << fObsLevel[i]/100. << "m";
252    *fLog << endl;
253
254    *fLog << "MagneticField:  X/Z=(" << fMagneticFieldX << "/";
255    *fLog << fMagneticFieldZ << ")" << UTF8::kMu << "T  Az=" << fMagneticFieldAz*TMath::RadToDeg() << UTF8::kDeg << "  (magnetic North w.r.t. North)" << endl;
256
257    *fLog << "Spectrum:       Slope=" << fSlopeSpectrum << "  (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" <<  endl;
258    *fLog << "Wavelength:     " << fWavelengthMin << "nm - " << fWavelengthMax << "nm" << endl;
259
260    if (fImpactMax>0)
261        *fLog << "ImpactMax:      " << fImpactMax << "cm" << endl;
262    if (fViewConeOuterAngle>0)
263        *fLog << "ViewCone:       " << fViewConeInnerAngle << UTF8::kDeg << " - " << fViewConeOuterAngle << UTF8::kDeg << endl;
264
265    *fLog << "Zd/Az:          ";
266    if (fZdMax>=0 && fZdMin<360)
267    {
268        *fLog << fZdMin << UTF8::kDeg;
269        if (fZdMin==fZdMax)
270            *fLog << " (fixed)";
271        else
272            *fLog << "-" << fZdMax << UTF8::kDeg;
273        *fLog << " / " << fAzMin << UTF8::kDeg;
274        if (fAzMin==fAzMax)
275            *fLog << " (fixed)";
276        else
277            *fLog << "-" << fAzMax << UTF8::kDeg;
278        *fLog << "  w.r.t. magnetic North." << endl;
279    }
280
281    if (fZdMin>=360) // 4010.7
282        *fLog << "-trajectory-" << endl;
283
284
285    if (fImpactMax>0)
286        *fLog << "Max.sim.Impact: " << fImpactMax << "cm" << endl;
287
288    *fLog << "Options used:  ";
289    if (Has(kCerenkov))
290        *fLog << " CERENKOV";
291    if (Has(kIact))
292        *fLog << " IACT";
293    if (Has(kCeffic))
294        *fLog << " CEFFIC";
295    if (Has(kAtmext))
296        *fLog << " ATMEXT" << GetNumAtmosphericModel();
297    if (Has(kRefraction))
298        *fLog << " +Refraction";
299    if (Has(kVolumedet))
300        *fLog << " VOLUMEDET";
301    if (Has(kCurved))
302        *fLog << " CURVED";
303    if (Has(kSlant))
304        *fLog << " SLANT";
305    *fLog << " [" << hex << fCerenkovFlag << "]" << dec << endl;
306
307    if (HasLayers())
308    {
309        *fLog << "Atm.Layers:    ";
310        for (int i=0; i<5; i++)
311            *fLog << " " << fAtmosphericLayers[i];
312    }
313    *fLog << endl;
314    *fLog << "Atm.Coeff A:   ";
315    for (int i=0; i<5; i++)
316        *fLog << " " << fAtmosphericCoeffA[i];
317    *fLog << endl;
318    *fLog << "Atm.Coeff B:   ";
319    for (int i=0; i<5; i++)
320        *fLog << " " << fAtmosphericCoeffB[i];
321    *fLog << endl;
322    *fLog << "Atm.Coeff C:   ";
323    for (int i=0; i<5; i++)
324        *fLog << " " << fAtmosphericCoeffC[i];
325    *fLog << endl;
326
327
328}
329
Note: See TracBrowser for help on using the repository browser.