source: trunk/MagicSoft/Mars/mhist/MHCamEventRot.cc@ 8389

Last change on this file since 8389 was 5158, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 8.1 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, 4/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHCamEventRot
28//
29// Derotate MCamEvent before filling a 2D histogram.
30//
31// Axis titles
32// ===========
33//
34// 1) If no other title is given 'a.u.' or 'counts'is used.
35// 2) If the corresponding MBinning has no default title. This title
36// is used for the z-axis
37// 3) If the title of MHCamEvent is different from the default,
38// it is used as histogram title. You can use this to set the
39// axis title, too. For more information see TH1::SetTitle, eg.
40// SetTitle("MyHist;;;z[photons]");
41// The x- and y-axis title is ignored.
42//
43// For example:
44// MHCamEventRot myhist("Title;;;z[photons]");
45//
46//
47// To be done: * Set the sky position of the center of the display and
48// correct if the pointing position of the telescope is
49// different.
50// * Make it quadratic like MHCamera
51//
52//////////////////////////////////////////////////////////////////////////////
53#include "MHCamEventRot.h"
54
55#include <TF1.h>
56#include <TH2.h>
57#include <TGraph.h>
58#include <TStyle.h>
59#include <TCanvas.h>
60#include <TPaveText.h>
61#include <TStopwatch.h>
62
63#include "MGeomCam.h"
64#include "MSrcPosCam.h"
65#include "MHillasSrc.h"
66#include "MTime.h"
67#include "MObservatory.h"
68#include "MPointingPos.h"
69#include "MAstroCatalog.h"
70#include "MAstroSky2Local.h"
71#include "MStatusDisplay.h"
72#include "MMath.h"
73
74#include "MBinning.h"
75#include "MParList.h"
76#include "MCamEvent.h"
77
78#include "MLog.h"
79#include "MLogManip.h"
80
81ClassImp(MHCamEventRot);
82
83using namespace std;
84
85const TString MHCamEventRot::gsDefName = "MHCamEventRot";
86const TString MHCamEventRot::gsDefTitle = "Plot showing derotated MCamEvent data";
87
88// --------------------------------------------------------------------------
89//
90// Default Constructor
91//
92MHCamEventRot::MHCamEventRot(const char *name, const char *title)
93 : fTime(0), fPointPos(0), fObservatory(0), fType(0), fNameTime("MTime"),
94 fThreshold(0), fUseThreshold(kNoBound)
95{
96 //
97 // set the name and title of this object
98 //
99 fName = name ? name : gsDefName.Data();
100 fTitle = title ? title : gsDefTitle.Data();
101
102 fHist.SetDirectory(NULL);
103
104 fHist.SetName("Derot");
105 fHist.SetTitle("2D-plot of MCamEvents (derotated)");
106}
107
108// --------------------------------------------------------------------------
109//
110// Set binnings (takes BinningCamEvent) and prepare filling of the
111// histogram.
112//
113// Also search for MTime, MObservatory and MPointingPos
114//
115Bool_t MHCamEventRot::SetupFill(const MParList *plist)
116{
117 fGeom = (MGeomCam*)plist->FindObject("MGeomCam");
118 if (!fGeom)
119 {
120 *fLog << err << "MGeomCam not found... aborting." << endl;
121 return kFALSE;
122 }
123
124 const MBinning *bins = (MBinning*)plist->FindObject("BinningCamEvent");
125 if (!bins)
126 {
127 const Float_t r = fGeom->GetMaxRadius()*fGeom->GetConvMm2Deg();
128
129 MBinning b;
130 b.SetEdges(41, -r, r);
131 SetBinning(&fHist, &b, &b);
132 }
133 else
134 SetBinning(&fHist, bins, bins);
135
136 fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
137 if (!fPointPos)
138 *fLog << warn << "MPointingPos not found... no derotation." << endl;
139
140 fTime = (MTime*)plist->FindObject(AddSerialNumber(fNameTime));
141 if (!fTime)
142 *fLog << warn << "MTime not found... no derotation." << endl;
143
144 fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
145 if (!fObservatory)
146 *fLog << warn << "MObservatory not found... no derotation." << endl;
147
148 // FIXME: Because the pointing position could change we must check
149 // for the current pointing position and add a offset in the
150 // Fill function!
151 if (fPointPos)
152 {
153 fRa = fPointPos->GetRa();
154 fDec = fPointPos->GetDec();
155 }
156
157 if (bins && bins->HasTitle())
158 fHist.SetZTitle(bins->GetTitle());
159 else
160 if (fTitle!=gsDefTitle)
161 {
162 fHist.SetTitle(fTitle);
163 fHist.SetXTitle("x [\\circ]");
164 fHist.SetYTitle("y [\\circ]");
165 }
166 else
167 if (!fTitle.Contains(";"))
168 fHist.SetZTitle(fUseThreshold==kNoBound?"a.u.":"Counts");
169
170 return kTRUE;
171}
172
173// --------------------------------------------------------------------------
174//
175// Fill the histogram. For details see the code or the class description
176//
177Bool_t MHCamEventRot::Fill(const MParContainer *par, const Stat_t w)
178{
179 const MCamEvent *evt = dynamic_cast<const MCamEvent*>(par);
180 if (!evt)
181 {
182 *fLog << err << "MHCamEventRot::Fill: No container specified!" << endl;
183 return kFALSE;
184 }
185
186 // Get max radius...
187 // const Double_t maxr = fGeom->GetMaxRadius()*fGeom->GetConvMm2Deg();
188
189 // Get camera rotation angle
190 Double_t rho = 0;
191 if (fTime && fObservatory && fPointPos)
192 rho = fPointPos->RotationAngle(*fObservatory, *fTime);
193 //if (fPointPos)
194 // rho = fPointPos->RotationAngle(*fObservatory);
195
196 // Get number of bins and bin-centers
197 const Int_t nx = fHist.GetNbinsX();
198 const Int_t ny = fHist.GetNbinsY();
199 Axis_t cx[nx];
200 Axis_t cy[ny];
201 fHist.GetXaxis()->GetCenter(cx);
202 fHist.GetYaxis()->GetCenter(cy);
203
204 for (int ix=0; ix<nx; ix++)
205 {
206 for (int iy=0; iy<ny; iy++)
207 {
208 // check distance... to get a circle plot
209 //if (TMath::Hypot(cx[ix], cy[iy])>maxr)
210 // continue;
211
212 // rotate center of bin
213 TVector2 v(cx[ix], cy[iy]);
214 if (rho!=0)
215 v=v.Rotate(rho);
216
217 // FIXME: slow! Create Lookup table instead!
218 const Int_t idx = fGeom->GetPixelIdxDeg(v);
219
220 Double_t val;
221 if (idx<0 || !evt->GetPixelContent(val, idx, *fGeom, fType))
222 continue;
223
224 // Fill histogram
225 if (fUseThreshold!=kNoBound)
226 {
227 if (val>fThreshold && fUseThreshold==kIsLowerBound ||
228 val<fThreshold && fUseThreshold==kIsUpperBound)
229 fHist.Fill(cx[ix], cy[iy]);
230 }
231 else
232 fHist.Fill(cx[ix], cy[iy], val);
233 }
234 }
235
236 return kTRUE;
237}
238
239// --------------------------------------------------------------------------
240//
241// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude
242// is set to 9, while the fov is equal to the current fov of the false
243// source plot.
244//
245TObject *MHCamEventRot::GetCatalog()
246{
247 const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
248
249 // Create catalog...
250 MAstroCatalog stars;
251 stars.SetLimMag(9);
252 stars.SetGuiActive(kFALSE);
253 stars.SetRadiusFOV(maxr);
254 stars.SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad());
255 stars.ReadBSC("bsc5.dat");
256
257 TObject *o = (MAstroCatalog*)stars.Clone();
258 o->SetBit(kCanDelete);
259
260 return o;
261}
262
263// --------------------------------------------------------------------------
264//
265// Draw the histogram
266//
267void MHCamEventRot::Draw(Option_t *opt)
268{
269 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
270 pad->SetBorderMode(0);
271
272 AppendPad("");
273
274 // draw the 2D histogram Sigmabar versus Theta
275 fHist.Draw(opt);
276 TObject *catalog = GetCatalog();
277 catalog->Draw("mirror same");
278}
Note: See TracBrowser for help on using the repository browser.