source: branches/MarsMoreSimulationTruth/mcorsika/MCorsikaEvtHeader.cc

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