source: trunk/Mars/mcorsika/MCorsikaEvtHeader.cc@ 20115

Last change on this file since 20115 was 19406, checked in by tbretz, 6 years ago
The version number had to be changed as the data member description changed (I think that is new...)
File size: 8.0 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 11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: Software Development, 2000-2010
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MCorsikaEvtHeader
28//
29// Class Version 2:
30// ----------------
31// - UInt_t fParticleID
32//
33//
34// Class Version 3:
35// ----------------
36// - UInt_t fNumReuse
37//
38//
39// Class Version 4:
40// ----------------
41// - Variable description changed
42//
43//
44/////////////////////////////////////////////////////////////////////////////
45#include "MCorsikaEvtHeader.h"
46#include "MCorsikaFormat.h"
47
48#include <iomanip>
49#include <fstream>
50
51#include "fits.h"
52
53#include <TMath.h>
54
55#include "MLog.h"
56#include "MLogManip.h"
57
58//#include "MMcEvt.hxx"
59
60ClassImp(MCorsikaEvtHeader);
61
62using namespace std;
63
64// --------------------------------------------------------------------------
65//
66// Default constructor. Create the array to store the data.
67//
68MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title)
69 : fNumReuse((UInt_t)-1), fX(0), fY(0)
70{
71 fName = name ? name : "MCorsikaEvtHeader";
72 fTitle = title ? title : "Raw Event Header Information";
73}
74
75// --------------------------------------------------------------------------
76//
77// Return Impact (distance of ground incident point from telescope axis)
78//
79// Distance d between a point q and a line u (via p):
80// d = | ( vec(q) - vec(p) ) x vec(u) | / | vec(u) |
81//
82// If p = 0
83//
84// ==> d = | vec(q) x vec(u) | / | vec(u) |
85// w := q = (x/y/0)
86// v := u = (Alt/Az)
87//
88// For Alt/Az = Theta/Phi:
89//
90// x = r sin(theta) cos(phi)
91// y = r sin(theta) sin(phi)
92// z = r cos(theta)
93//
94// ( -q3*u2 ) ( -cos(theta)*y )
95// vec(q) x vec(u) = ( q3*u1 ) = ( cos(theta)*x )
96// ( q1*u2-q2*u1 ) ( sin(theta) cos(phi) * y - sin(theta) sin(phi) * x )
97//
98// ==> d = sqrt( cos(theta)^2 (x^2 + y^2 ) +
99// sin(theta)^2 ( cos(phi) y + sin(phi) x)^2 )
100//
101Double_t MCorsikaEvtHeader::GetImpact() const
102{
103 const Double_t c = TMath::Cos(fZd);
104 const Double_t s = TMath::Sin(fZd);
105
106 const Double_t p = TMath::Cos(fAz)*fX - TMath::Sin(fAz)*fY;
107
108 return TMath::Sqrt(c*c*(fX*fX + fY*fY) + s*s* p*p);
109}
110
111// --------------------------------------------------------------------------
112//
113// This member function prints all Data of one Event to *fLog.
114//
115void MCorsikaEvtHeader::Print(Option_t *o) const
116{
117 *fLog << all;
118 *fLog << "Event Number: " << dec << fEvtNumber << " (reused=" << fNumReuse << ")" << endl;
119// *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl;
120 *fLog << "Energy: " << fTotalEnergy << "GeV" << endl;
121 *fLog << "Starting Altitude: " << fStartAltitude << "g/cm" << UTF8::kSquare << endl;
122 *fLog << "Number of 1st Target: " << fFirstTargetNum << endl,
123 *fLog << "Height of 1st Interaction: " << fFirstInteractionHeight/100. << "m" << endl;
124 *fLog << "Momentum X/Y/Z (GeV/c): " << fMomentumX << "/" << fMomentumY << "/" << fMomentumZ << endl;
125 *fLog << "Zenith/Azimuth Angle: " << fZd*TMath::RadToDeg() << UTF8::kDeg << "/" << fAz*TMath::RadToDeg() << UTF8::kDeg << endl;
126 *fLog << "Impact X/Y: " << fX/100. << "m/" << fY/100. << "m (r=" << TMath::Hypot(fX, fY)/100. << "m)" << endl;
127 *fLog << "Weighted Num Photons: " << fWeightedNumPhotons << endl;
128 *fLog << endl;
129}
130
131// --------------------------------------------------------------------------
132//
133// get the EVTH information from the input parameter f
134//
135Int_t MCorsikaEvtHeader::ReadEvt(Float_t * f)
136{
137
138 fEvtNumber = TMath::Nint(f[0]);
139// fParticleID = TMath::Nint(f[1]);
140
141 fTotalEnergy = f[2];
142 fStartAltitude = f[3];
143 fFirstTargetNum = f[4];
144 fFirstInteractionHeight = f[5]; // FIXME: Could be negative
145
146 // Pointing opposite particle direction
147 // (x=north, y=west, z=upwards)
148 //fMomentumX = f[6];
149 //fMomentumY = f[7];
150 //fMomentumZ = f[8];
151 // Pointing along particle direction
152 // (x=east, y=north, z=upwards)
153 fMomentumX = f[7];
154 fMomentumY = -f[6];
155 fMomentumZ = -f[8]; // Does this minus make sense?!
156
157 fZd = f[9];
158 fAz = TMath::Pi()-f[10];
159
160 // f[96] // Number of reuse of event [max=20]
161 fTotReuse = TMath::Nint(f[96]);
162 if (fTotReuse > 20)
163 {
164 *fLog << err << "Number of reuse of shower is " << fTotReuse
165 << " But maximum implemented: 20" << endl;
166 return kERROR;
167 }
168
169 memcpy(fTempX, f +97, 20*sizeof(Float_t));
170 memcpy(fTempY, f+117, 20*sizeof(Float_t));
171
172
173 // FIXME: Check fNumReuse<20
174// fX = f[117 + fNumReuse];
175// fY = -f[ 97 + fNumReuse];
176
177 fWeightedNumPhotons = 0;
178
179 return kTRUE;
180}
181
182
183// this member function is for reading the event end block
184
185Int_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat, const uint32_t &blockLength)
186{
187 vector<Float_t> f(blockLength);
188 if (!fInFormat->Read(f.data(), blockLength*sizeof(Float_t)))
189 return kERROR;
190
191 const UInt_t evtnum = TMath::Nint(f[0]);
192 if (evtnum!=fEvtNumber)
193 {
194 *fLog << err << "ERROR - Mismatch in stream: Event number in EVTE (";
195 *fLog << evtnum << ") doesn't match EVTH (" << fEvtNumber << ")." << endl;
196 return kERROR;
197 }
198
199 // FIXME: What's the meaning of this if in case of reusing the event this number
200 // exists only once?
201 fWeightedNumPhotons = f[1];
202
203 return fInFormat->Eof() ? kERROR : kTRUE;
204}
205
206Bool_t MCorsikaEvtHeader::SetupFits(fits &fin)
207{
208 if (!fin.SetRefAddress("MCorsikaEvtHeader.fEvtNumber", fEvtNumber )) return kFALSE;
209 if (!fin.SetRefAddress("MCorsikaEvtHeader.fNumReuse", fNumReuse )) return kFALSE;
210 if (!fin.SetRefAddress("MCorsikaEvtHeader.fTotalEnergy", fTotalEnergy )) return kFALSE;
211 if (!fin.SetRefAddress("MCorsikaEvtHeader.fStartAltitude", fStartAltitude )) return kFALSE;
212 if (!fin.SetRefAddress("MCorsikaEvtHeader.fFirstTargetNum", fFirstTargetNum )) return kFALSE;
213 if (!fin.SetRefAddress( "CorsikaEvtHeader.fFirstInteractionHeight", fFirstInteractionHeight )) return kFALSE;
214 if (!fin.SetRefAddress("MCorsikaEvtHeader.fMomentumX", fMomentumX )) return kFALSE;
215 if (!fin.SetRefAddress("MCorsikaEvtHeader.fMomentumY", fMomentumY )) return kFALSE;
216 if (!fin.SetRefAddress("MCorsikaEvtHeader.fMomentumZ", fMomentumZ )) return kFALSE;
217 if (!fin.SetRefAddress("MCorsikaEvtHeader.fZd", fZd )) return kFALSE;
218 if (!fin.SetRefAddress("MCorsikaEvtHeader.fAz", fAz )) return kFALSE;
219 if (!fin.SetRefAddress("MCorsikaEvtHeader.fX", fX )) return kFALSE;
220 if (!fin.SetRefAddress("MCorsikaEvtHeader.fY", fY )) return kFALSE;
221 if (!fin.SetRefAddress("MCorsikaEvtHeader.fWeightedNumPhotons", fWeightedNumPhotons )) return kFALSE;
222
223 return kTRUE;
224}
Note: See TracBrowser for help on using the repository browser.