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):
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2020
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | /////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MMcEvt
|
---|
28 | //
|
---|
29 | // This class handles and contains the MonteCarlo information
|
---|
30 | // with which the events have been generated
|
---|
31 | // This information exists for each event.
|
---|
32 | //
|
---|
33 | // Note: All the azimuth Phi angles in this and other MC classes follow
|
---|
34 | // the convention in the Corsika program (see Corsika manual). There, phi
|
---|
35 | // is the azimuth of the momentum vector of particles, and is measured
|
---|
36 | // from the north direction, anticlockwise (i.e, west is phi=90 degrees).
|
---|
37 | // When it refers to the telescope orientation, it is the azimuth of a
|
---|
38 | // vector along the telescope axis, going from the camera to the mirror.
|
---|
39 | // So, fTelescopeTheta=90, fTelescopePhi = 0 means the telescope is
|
---|
40 | // pointing horizontally towards South. For an explanation, see also
|
---|
41 | // TDAS 02-11.
|
---|
42 | //
|
---|
43 | // Version 4:
|
---|
44 | // - Added member fFadcTimeJitter
|
---|
45 | //
|
---|
46 | // Version 5:
|
---|
47 | // - removed fPartId, fEnergy, fImpact, fTelescopeTheta, fTelescopePhi
|
---|
48 | // - derives now from MMcEvtBasic which contains all these values
|
---|
49 | // - moved ParticleId_t to base class MMcEvtBasic
|
---|
50 | //
|
---|
51 | // Version 6:
|
---|
52 | // - added fEventReuse
|
---|
53 | //
|
---|
54 | // Version 7:
|
---|
55 | // - moved fPhi to MMcEvtBasic
|
---|
56 | // - moved fTheta to MmcEvtBasic
|
---|
57 | //
|
---|
58 | // Version 8:
|
---|
59 | // - added fRunNumber
|
---|
60 | //
|
---|
61 | // Version 9:
|
---|
62 | // - moved CoreX to MMcEvtBasic
|
---|
63 | // - moved CoreY to MMcEvtBasic
|
---|
64 | //
|
---|
65 | /////////////////////////////////////////////////////////////////////////////
|
---|
66 | #include "MMcEvt.hxx"
|
---|
67 |
|
---|
68 | #include "fits.h"
|
---|
69 |
|
---|
70 | #include <TMath.h>
|
---|
71 |
|
---|
72 | #include "MLog.h"
|
---|
73 | #include "MLogManip.h"
|
---|
74 |
|
---|
75 | #include "MString.h"
|
---|
76 |
|
---|
77 | ClassImp(MMcEvt);
|
---|
78 |
|
---|
79 | using namespace std;
|
---|
80 |
|
---|
81 | // --------------------------------------------------------------------------
|
---|
82 | //
|
---|
83 | // Default constructor. Calls Clear()
|
---|
84 | //
|
---|
85 | MMcEvt::MMcEvt()
|
---|
86 | {
|
---|
87 | fName = "MMcEvt";
|
---|
88 | fTitle = "Event info from Monte Carlo";
|
---|
89 |
|
---|
90 | Clear();
|
---|
91 | }
|
---|
92 |
|
---|
93 | // --------------------------------------------------------------------------
|
---|
94 | //
|
---|
95 | // Constructor. Use this to set all data members
|
---|
96 | //
|
---|
97 | // THIS FUNCTION IS FOR THE SIMULATION OLNY.
|
---|
98 | // DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS.
|
---|
99 | //
|
---|
100 | MMcEvt::MMcEvt(UInt_t fEvtNum, ParticleId_t usPId, Float_t fEner,
|
---|
101 | Float_t fThi0, Float_t fFirTar, Float_t fzFirInt,
|
---|
102 | Float_t fThet, Float_t fPhii, Float_t fCorD,
|
---|
103 | Float_t fCorX, Float_t fCorY, Float_t fImpa,
|
---|
104 | Float_t fTPhii, Float_t fTThet, Float_t fTFirst,
|
---|
105 | Float_t fTLast, Float_t fL_Nmax, Float_t fL_t0,
|
---|
106 | Float_t fL_tmax, Float_t fL_a, Float_t fL_b,
|
---|
107 | Float_t fL_c, Float_t fL_chi2, UInt_t uiPin,
|
---|
108 | UInt_t uiPat, UInt_t uiPre, UInt_t uiPco,
|
---|
109 | UInt_t uiPelS, UInt_t uiPelC, Float_t elec,
|
---|
110 | Float_t muon, Float_t other, Float_t fadc_jitter,
|
---|
111 | Int_t reuse) : fRunNumber(0)
|
---|
112 | {
|
---|
113 | fName = "MMcEvt";
|
---|
114 | fTitle = "Event info from Monte Carlo";
|
---|
115 |
|
---|
116 | Fill(fEvtNum, usPId, fEner, fThi0, fFirTar, fzFirInt, fThet,
|
---|
117 | fPhii, fCorD, fCorX, fCorY, fImpa, fTPhii, fTThet, fTFirst,
|
---|
118 | fTLast, fL_Nmax, fL_t0, fL_tmax, fL_a, fL_b, fL_c, fL_chi2,
|
---|
119 | uiPin, uiPat, uiPre, uiPco, uiPelS, uiPelC, elec, muon, other,
|
---|
120 | fadc_jitter, reuse);
|
---|
121 | }
|
---|
122 |
|
---|
123 | // --------------------------------------------------------------------------
|
---|
124 | //
|
---|
125 | // reset all values to values as nonsense as possible
|
---|
126 | //
|
---|
127 | void MMcEvt::Clear(Option_t *opt)
|
---|
128 | {
|
---|
129 | MMcEvtBasic::Clear();
|
---|
130 |
|
---|
131 | fRunNumber = 0;
|
---|
132 |
|
---|
133 | fCoreD = 0;
|
---|
134 |
|
---|
135 | fPhotIni = 0;
|
---|
136 | fPassPhotAtm = 0;
|
---|
137 | fPassPhotRef = 0;
|
---|
138 | fPassPhotCone = 0;
|
---|
139 | fPhotElfromShower = 0;
|
---|
140 | fPhotElinCamera = 0;
|
---|
141 |
|
---|
142 | fElecCphFraction=0;
|
---|
143 | fMuonCphFraction=0;
|
---|
144 | fOtherCphFraction=0;
|
---|
145 |
|
---|
146 | fEventReuse = 0;
|
---|
147 | }
|
---|
148 |
|
---|
149 | // --------------------------------------------------------------------------
|
---|
150 | //
|
---|
151 | // Use this to set all data members
|
---|
152 | //
|
---|
153 | // THIS FUNCTION IS FOR THE SIMULATION OLNY.
|
---|
154 | // DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS.
|
---|
155 | //
|
---|
156 | void MMcEvt::Fill( UInt_t fEvtNum, ParticleId_t usPId, Float_t fEner,
|
---|
157 | Float_t fThi0, Float_t fFirTar, Float_t fzFirInt,
|
---|
158 | Float_t fThet, Float_t fPhii, Float_t fCorD,
|
---|
159 | Float_t fCorX, Float_t fCorY, Float_t fImpa,
|
---|
160 | Float_t fTPhii, Float_t fTThet, Float_t fTFirst,
|
---|
161 | Float_t fTLast, Float_t fL_Nmax, Float_t fL_t0,
|
---|
162 | Float_t fL_tmax, Float_t fL_a, Float_t fL_b,
|
---|
163 | Float_t fL_c, Float_t fL_chi2, UInt_t uiPin,
|
---|
164 | UInt_t uiPat, UInt_t uiPre, UInt_t uiPco,
|
---|
165 | UInt_t uiPelS, UInt_t uiPelC, Float_t elec,
|
---|
166 | Float_t muon, Float_t other, Float_t fadc_jitter,
|
---|
167 | Int_t reuse)
|
---|
168 | {
|
---|
169 | //
|
---|
170 | // All datamembers are filled with the correspondin parameters.
|
---|
171 | //
|
---|
172 | // Don't use this memberfunction in analysis
|
---|
173 | //
|
---|
174 | fEvtNumber = fEvtNum;
|
---|
175 | fPartId = usPId ;
|
---|
176 | fEnergy = fEner ;
|
---|
177 | fThick0 = fThi0;
|
---|
178 | fFirstTarget = fFirTar;
|
---|
179 | fZFirstInteraction = fzFirInt;
|
---|
180 |
|
---|
181 | fTheta = fThet ;
|
---|
182 | fPhi = fPhii ;
|
---|
183 |
|
---|
184 | fCoreD = fCorD ;
|
---|
185 | fCoreX = fCorX ;
|
---|
186 | fCoreY = fCorY ;
|
---|
187 | fImpact = fImpa ;
|
---|
188 |
|
---|
189 | fTelescopePhi = fTPhii;
|
---|
190 | fTelescopeTheta = fTThet;
|
---|
191 | fTimeFirst = fTFirst;
|
---|
192 | fTimeLast = fTLast;
|
---|
193 | fLongiNmax = fL_Nmax;
|
---|
194 | fLongit0 = fL_t0;
|
---|
195 | fLongitmax = fL_tmax;
|
---|
196 | fLongia = fL_a;
|
---|
197 | fLongib = fL_b;
|
---|
198 | fLongic = fL_c;
|
---|
199 | fLongichi2 = fL_chi2;
|
---|
200 |
|
---|
201 | fPhotIni = uiPin ;
|
---|
202 | fPassPhotAtm = fPhotIni-uiPat ;
|
---|
203 | fPassPhotRef = fPassPhotAtm-uiPre ;
|
---|
204 | fPassPhotCone = uiPco ;
|
---|
205 | fPhotElfromShower = uiPelS ;
|
---|
206 | fPhotElinCamera = uiPelC ;
|
---|
207 |
|
---|
208 | fElecCphFraction=elec;
|
---|
209 | fMuonCphFraction=muon;
|
---|
210 | fOtherCphFraction=other;
|
---|
211 |
|
---|
212 | fFadcTimeJitter = fadc_jitter;
|
---|
213 |
|
---|
214 | fEventReuse = reuse;
|
---|
215 | }
|
---|
216 |
|
---|
217 | // --------------------------------------------------------------------------
|
---|
218 | //
|
---|
219 | // Print the contents of the container.
|
---|
220 | //
|
---|
221 | // if you specify an option only the requested data members are printed:
|
---|
222 | // allowed options are:
|
---|
223 | // id, energy, impact
|
---|
224 | //
|
---|
225 | void MMcEvt::Print(Option_t *opt) const
|
---|
226 | {
|
---|
227 | MMcEvtBasic::Print(opt);
|
---|
228 |
|
---|
229 | TString str(opt);
|
---|
230 | if (str.IsNull())
|
---|
231 | {
|
---|
232 | *fLog << " ID: Run=" << fRunNumber << ", Evt=" << fEvtNumber << endl;
|
---|
233 | *fLog << " Photoelectrons: " << fPhotElfromShower << endl;
|
---|
234 | }
|
---|
235 | }
|
---|
236 |
|
---|
237 | // --------------------------------------------------------------------------
|
---|
238 | //
|
---|
239 | // Return a proper description of the monte carlo event
|
---|
240 | //
|
---|
241 | TString MMcEvt::GetDescription(const TString &s) const
|
---|
242 | {
|
---|
243 | TString txt("#splitline{");
|
---|
244 |
|
---|
245 | txt += GetParticleName();
|
---|
246 | txt += " ";
|
---|
247 | txt += s;
|
---|
248 | txt += "}{ E=";
|
---|
249 | txt += GetEnergyStr();
|
---|
250 | txt += " r=";
|
---|
251 | txt += TMath::Nint(GetImpact()/100);
|
---|
252 | txt += "m Zd=";
|
---|
253 | txt += MString::Format("%.1f", GetTelescopeTheta()*TMath::RadToDeg());
|
---|
254 | txt += "\\circ ";
|
---|
255 | if (GetPhotElfromShower()>=10000)
|
---|
256 | txt += MString::Format("%dk", TMath::Nint(GetPhotElfromShower()/1000.));
|
---|
257 | else
|
---|
258 | if (GetPhotElfromShower()>=1000)
|
---|
259 | txt += MString::Format("%.1fk", GetPhotElfromShower()/1000.);
|
---|
260 | else
|
---|
261 | txt += GetPhotElfromShower();
|
---|
262 | txt += "PhEl}";
|
---|
263 |
|
---|
264 | return txt;
|
---|
265 | }
|
---|
266 |
|
---|
267 | Bool_t MMcEvt::SetupFits(fits &fin)
|
---|
268 | {
|
---|
269 | if (fName!="MMcEvt")
|
---|
270 | {
|
---|
271 | *fLog << err << "ERROR - SetupFits only supported if name equals MMcEvt." << endl;
|
---|
272 | return kFALSE;
|
---|
273 | }
|
---|
274 |
|
---|
275 | if (!fin.SetRefAddress("MMcEvt.fEvtNumber", fEvtNumber )) return kFALSE;
|
---|
276 | if (!fin.SetRefAddress("MMcEvt.fThick0", fThick0 )) return kFALSE;
|
---|
277 | if (!fin.SetRefAddress("MMcEvt.fFirstTarget", fFirstTarget )) return kFALSE;
|
---|
278 | if (!fin.SetRefAddress("MMcEvt.fZFirstInteraction", fZFirstInteraction)) return kFALSE;
|
---|
279 | if (!fin.SetRefAddress("MMcEvt.fCoreD", fCoreD )) return kFALSE;
|
---|
280 | if (!fin.SetRefAddress("MMcEvt.fTimeFirst", fTimeFirst )) return kFALSE;
|
---|
281 | if (!fin.SetRefAddress("MMcEvt.fTimeLast", fTimeLast )) return kFALSE;
|
---|
282 | if (!fin.SetRefAddress("MMcEvt.fLongiNmax", fLongiNmax )) return kFALSE;
|
---|
283 | if (!fin.SetRefAddress("MMcEvt.fLongit0", fLongit0 )) return kFALSE;
|
---|
284 | if (!fin.SetRefAddress("MMcEvt.fLongitmax", fLongitmax )) return kFALSE;
|
---|
285 | if (!fin.SetRefAddress("MMcEvt.fLongia", fLongia )) return kFALSE;
|
---|
286 | if (!fin.SetRefAddress("MMcEvt.fLongib", fLongib )) return kFALSE;
|
---|
287 | if (!fin.SetRefAddress("MMcEvt.fLongic", fLongic )) return kFALSE;
|
---|
288 | if (!fin.SetRefAddress("MMcEvt.fLongichi2", fLongichi2 )) return kFALSE;
|
---|
289 | if (!fin.SetRefAddress("MMcEvt.fPhotIni", fPhotIni )) return kFALSE;
|
---|
290 | if (!fin.SetRefAddress("MMcEvt.fPassPhotAtm", fPassPhotAtm )) return kFALSE;
|
---|
291 | if (!fin.SetRefAddress("MMcEvt.fPassPhotRef", fPassPhotRef )) return kFALSE;
|
---|
292 | if (!fin.SetRefAddress("MMcEvt.fPassPhotCone", fPassPhotCone )) return kFALSE;
|
---|
293 | if (!fin.SetRefAddress("MMcEvt.fPhotElfromShower", fPhotElfromShower )) return kFALSE;
|
---|
294 | if (!fin.SetRefAddress("MMcEvt.fPhotElinCamera", fPhotElinCamera )) return kFALSE;
|
---|
295 | if (!fin.SetRefAddress("MMcEvt.fElecCphFraction", fElecCphFraction )) return kFALSE;
|
---|
296 | if (!fin.SetRefAddress("MMcEvt.fMuonCphFraction", fMuonCphFraction )) return kFALSE;
|
---|
297 | if (!fin.SetRefAddress("MMcEvt.fOtherCphFraction", fOtherCphFraction )) return kFALSE;
|
---|
298 | if (!fin.SetRefAddress("MMcEvt.fFadcTimeJitter", fFadcTimeJitter )) return kFALSE;
|
---|
299 | if (!fin.SetRefAddress("MMcEvt.fEventReuse", fEventReuse )) return kFALSE;
|
---|
300 |
|
---|
301 | // Not in files written from MMcEvt < V8
|
---|
302 | if (fin.HasColumn("MMcEvt.fRunNumber"))
|
---|
303 | if (!fin.SetRefAddress("MMcEvt.fRunNumber", fRunNumber))
|
---|
304 | return kFALSE;
|
---|
305 |
|
---|
306 | // Written to MmcEvtBasic for MMcEvt > V8
|
---|
307 | if (fin.HasColumn("MMcEvt.fCoreX"))
|
---|
308 | if (!fin.SetRefAddress("MMcEvt.fCoreX", fCoreX))
|
---|
309 | return kFALSE;
|
---|
310 |
|
---|
311 | if (fin.HasColumn("MMcEvt.fCoreY"))
|
---|
312 | if (!fin.SetRefAddress("MMcEvt.fCoreY", fCoreY))
|
---|
313 | return kFALSE;
|
---|
314 |
|
---|
315 | return MMcEvtBasic::SetupFits(fin);
|
---|
316 | }
|
---|