1 | /* ======================================================================== *\
|
---|
2 | ! $Name: not supported by cvs2svn $:$Id: MHCamera.cc,v 1.98 2007-02-01 11:24:18 tbretz Exp $
|
---|
3 | ! --------------------------------------------------------------------------
|
---|
4 | !
|
---|
5 | ! *
|
---|
6 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
7 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
8 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
9 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
10 | ! *
|
---|
11 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
12 | ! * documentation for any purpose is hereby granted without fee,
|
---|
13 | ! * provided that the above copyright notice appear in all copies and
|
---|
14 | ! * that both that copyright notice and this permission notice appear
|
---|
15 | ! * in supporting documentation. It is provided "as is" without express
|
---|
16 | ! * or implied warranty.
|
---|
17 | ! *
|
---|
18 | !
|
---|
19 | !
|
---|
20 | ! Author(s): Thomas Bretz, 05/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
21 | ! Author(s): Harald Kornmayer, 1/2001
|
---|
22 | ! Author(s): Markus Gaug, 03/2004 <mailto:markus@ifae.es>
|
---|
23 | !
|
---|
24 | ! Copyright: MAGIC Software Development, 2000-2004
|
---|
25 | !
|
---|
26 | !
|
---|
27 | \* ======================================================================== */
|
---|
28 |
|
---|
29 | /////////////////////////////////////////////////////////////////////////////
|
---|
30 | //
|
---|
31 | // MHCamera
|
---|
32 | //
|
---|
33 | // Camera Display, based on a TH1D. Pleas be carefull using the
|
---|
34 | // underlaying TH1D.
|
---|
35 | //
|
---|
36 | // To change the scale to a logarithmic scale SetLogy() of the Pad.
|
---|
37 | //
|
---|
38 | // You can correct for the abberation. Assuming that the distance
|
---|
39 | // between the mean position of the light distribution and the position
|
---|
40 | // of a perfect reflection on a perfect mirror in the distance r on
|
---|
41 | // the camera plane is dr it is d = a*dr while a is the abberation
|
---|
42 | // constant (for the MAGIC mirror it is roughly 0.0713). You can
|
---|
43 | // set this constant by calling SetAbberation(a) which results in a
|
---|
44 | // 'corrected' display (all outer pixels are shifted towards the center
|
---|
45 | // of the camera to correct for this abberation)
|
---|
46 | //
|
---|
47 | // Be carefull: Entries in this context means Entries/bin or Events
|
---|
48 | //
|
---|
49 | // FIXME? Maybe MHCamera can take the fLog object from MGeomCam?
|
---|
50 | //
|
---|
51 | ////////////////////////////////////////////////////////////////////////////
|
---|
52 | #include "MHCamera.h"
|
---|
53 |
|
---|
54 | #include <fstream>
|
---|
55 | #include <iostream>
|
---|
56 |
|
---|
57 | #include <TBox.h>
|
---|
58 | #include <TArrow.h>
|
---|
59 | #include <TLatex.h>
|
---|
60 | #include <TStyle.h>
|
---|
61 | #include <TCanvas.h>
|
---|
62 | #include <TArrayF.h>
|
---|
63 | #include <TRandom.h>
|
---|
64 | #include <TPaveText.h>
|
---|
65 | #include <TPaveStats.h>
|
---|
66 | #include <TClonesArray.h>
|
---|
67 | #include <THistPainter.h>
|
---|
68 | #include <THLimitsFinder.h>
|
---|
69 | #include <TProfile.h>
|
---|
70 | #include <TH1.h>
|
---|
71 | #include <TF1.h>
|
---|
72 | #include <TCanvas.h>
|
---|
73 | #include <TLegend.h>
|
---|
74 |
|
---|
75 | #include "MLog.h"
|
---|
76 | #include "MLogManip.h"
|
---|
77 |
|
---|
78 | #include "MH.h"
|
---|
79 | #include "MBinning.h"
|
---|
80 | #include "MHexagon.h"
|
---|
81 |
|
---|
82 | #include "MGeomPix.h"
|
---|
83 | #include "MGeomCam.h"
|
---|
84 |
|
---|
85 | #include "MCamEvent.h"
|
---|
86 |
|
---|
87 | #include "MArrayD.h"
|
---|
88 | #include "MMath.h" // MMath::GaussProb
|
---|
89 |
|
---|
90 | #define kItemsLegend 48 // see SetPalette(1,0)
|
---|
91 |
|
---|
92 | ClassImp(MHCamera);
|
---|
93 |
|
---|
94 | using namespace std;
|
---|
95 |
|
---|
96 | void MHCamera::Init()
|
---|
97 | {
|
---|
98 | Sumw2();
|
---|
99 |
|
---|
100 | UseCurrentStyle();
|
---|
101 |
|
---|
102 | SetDirectory(NULL);
|
---|
103 |
|
---|
104 | SetLineColor(kGreen);
|
---|
105 | SetMarkerStyle(kFullDotMedium);
|
---|
106 | SetXTitle("Pixel Index");
|
---|
107 |
|
---|
108 | fNotify = new TList;
|
---|
109 | fNotify->SetBit(kMustCleanup);
|
---|
110 | gROOT->GetListOfCleanups()->Add(fNotify);
|
---|
111 |
|
---|
112 | TVirtualPad *save = gPad;
|
---|
113 | gPad = 0;
|
---|
114 | /*
|
---|
115 | #if ROOT_VERSION_CODE < ROOT_VERSION(3,01,06)
|
---|
116 | SetPalette(1, 0);
|
---|
117 | #endif
|
---|
118 | */
|
---|
119 | /*
|
---|
120 | #if ROOT_VERSION_CODE < ROOT_VERSION(4,04,00)
|
---|
121 | SetPrettyPalette();
|
---|
122 | #elese
|
---|
123 | // WORAROUND - FIXME: Calling it many times becomes slower and slower
|
---|
124 | SetInvDeepBlueSeaPalette();
|
---|
125 | #endif
|
---|
126 | */
|
---|
127 | gPad = save;
|
---|
128 | }
|
---|
129 |
|
---|
130 | // ------------------------------------------------------------------------
|
---|
131 | //
|
---|
132 | // Default Constructor. To be used by the root system ONLY.
|
---|
133 | //
|
---|
134 | MHCamera::MHCamera() : TH1D(), fGeomCam(NULL), fAbberation(0)
|
---|
135 | {
|
---|
136 | Init();
|
---|
137 | }
|
---|
138 |
|
---|
139 | // ------------------------------------------------------------------------
|
---|
140 | //
|
---|
141 | // Constructor. Makes a clone of MGeomCam. Removed the TH1D from the
|
---|
142 | // current directory. Calls Sumw2(). Set the histogram line color
|
---|
143 | // (for error bars) to Green and the marker style to kFullDotMedium.
|
---|
144 | //
|
---|
145 | MHCamera::MHCamera(const MGeomCam &geom, const char *name, const char *title)
|
---|
146 | : fGeomCam(NULL), fAbberation(0)
|
---|
147 | {
|
---|
148 | //fGeomCam = (MGeomCam*)geom.Clone();
|
---|
149 | SetGeometry(geom, name, title);
|
---|
150 | Init();
|
---|
151 |
|
---|
152 | //
|
---|
153 | // root 3.02
|
---|
154 | // * base/inc/TObject.h:
|
---|
155 | // register BIT(8) as kNoContextMenu. If an object has this bit set it will
|
---|
156 | // not get an automatic context menu when clicked with the right mouse button.
|
---|
157 | }
|
---|
158 |
|
---|
159 | void MHCamera::SetGeometry(const MGeomCam &geom, const char *name, const char *title)
|
---|
160 | {
|
---|
161 | SetNameTitle(name, title);
|
---|
162 |
|
---|
163 | TAxis &x = *GetXaxis();
|
---|
164 |
|
---|
165 | SetBins(geom.GetNumPixels(), 0, 1);
|
---|
166 | x.Set(geom.GetNumPixels(), -0.5, geom.GetNumPixels()-0.5);
|
---|
167 |
|
---|
168 | //SetBins(geom.GetNumPixels(), -0.5, geom.GetNumPixels()-0.5);
|
---|
169 | //Rebuild();
|
---|
170 |
|
---|
171 | if (fGeomCam)
|
---|
172 | delete fGeomCam;
|
---|
173 | fGeomCam = (MGeomCam*)geom.Clone();
|
---|
174 |
|
---|
175 | fUsed.Set(geom.GetNumPixels());
|
---|
176 | for (Int_t i=0; i<fNcells-2; i++)
|
---|
177 | ResetUsed(i);
|
---|
178 |
|
---|
179 | fBinEntries.Set(geom.GetNumPixels()+2);
|
---|
180 | fBinEntries.Reset();
|
---|
181 | }
|
---|
182 |
|
---|
183 | // ------------------------------------------------------------------------
|
---|
184 | //
|
---|
185 | // Destructor. Deletes the cloned fGeomCam and the notification list.
|
---|
186 | //
|
---|
187 | MHCamera::~MHCamera()
|
---|
188 | {
|
---|
189 | if (fGeomCam)
|
---|
190 | delete fGeomCam;
|
---|
191 | if (fNotify)
|
---|
192 | delete fNotify;
|
---|
193 | }
|
---|
194 |
|
---|
195 | // ------------------------------------------------------------------------
|
---|
196 | //
|
---|
197 | // Return kTRUE for sector<0. Otherwise return kTRUE only if the specified
|
---|
198 | // sector idx matches the sector of the pixel with index idx.
|
---|
199 | //
|
---|
200 | Bool_t MHCamera::MatchSector(Int_t idx, const TArrayI §or, const TArrayI &aidx) const
|
---|
201 | {
|
---|
202 | const MGeomPix &pix = (*fGeomCam)[idx];
|
---|
203 | return FindVal(sector, pix.GetSector()) && FindVal(aidx, pix.GetAidx());
|
---|
204 | }
|
---|
205 |
|
---|
206 | // ------------------------------------------------------------------------
|
---|
207 | //
|
---|
208 | // Taken from TH1D::Fill(). Uses the argument directly as bin index.
|
---|
209 | // Doesn't increment the number of entries.
|
---|
210 | //
|
---|
211 | // -*-*-*-*-*-*-*-*Increment bin with abscissa X by 1*-*-*-*-*-*-*-*-*-*-*
|
---|
212 | // ==================================
|
---|
213 | //
|
---|
214 | // if x is less than the low-edge of the first bin, the Underflow bin is incremented
|
---|
215 | // if x is greater than the upper edge of last bin, the Overflow bin is incremented
|
---|
216 | //
|
---|
217 | // If the storage of the sum of squares of weights has been triggered,
|
---|
218 | // via the function Sumw2, then the sum of the squares of weights is incremented
|
---|
219 | // by 1 in the bin corresponding to x.
|
---|
220 | //
|
---|
221 | // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
|
---|
222 | Int_t MHCamera::Fill(Axis_t x)
|
---|
223 | {
|
---|
224 | #if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
|
---|
225 | if (fBuffer) return BufferFill(x,1);
|
---|
226 | #endif
|
---|
227 | const Int_t bin = (Int_t)x+1;
|
---|
228 | AddBinContent(bin);
|
---|
229 | fBinEntries[bin]++;
|
---|
230 | if (fSumw2.fN)
|
---|
231 | fSumw2.fArray[bin]++;
|
---|
232 |
|
---|
233 | if (bin<=0 || bin>fNcells-2)
|
---|
234 | return -1;
|
---|
235 |
|
---|
236 | fTsumw++;
|
---|
237 | fTsumw2++;
|
---|
238 | fTsumwx += x;
|
---|
239 | fTsumwx2 += x*x;
|
---|
240 | return bin;
|
---|
241 | }
|
---|
242 |
|
---|
243 | // ------------------------------------------------------------------------
|
---|
244 | //
|
---|
245 | // Taken from TH1D::Fill(). Uses the argument directly as bin index.
|
---|
246 | // Doesn't increment the number of entries.
|
---|
247 | //
|
---|
248 | // -*-*-*-*-*-*Increment bin with abscissa X with a weight w*-*-*-*-*-*-*-*
|
---|
249 | // =============================================
|
---|
250 | //
|
---|
251 | // if x is less than the low-edge of the first bin, the Underflow bin is incremented
|
---|
252 | // if x is greater than the upper edge of last bin, the Overflow bin is incremented
|
---|
253 | //
|
---|
254 | // If the storage of the sum of squares of weights has been triggered,
|
---|
255 | // via the function Sumw2, then the sum of the squares of weights is incremented
|
---|
256 | // by w^2 in the bin corresponding to x.
|
---|
257 | //
|
---|
258 | // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
|
---|
259 | Int_t MHCamera::Fill(Axis_t x, Stat_t w)
|
---|
260 | {
|
---|
261 | #if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
|
---|
262 | if (fBuffer) return BufferFill(x,w);
|
---|
263 | #endif
|
---|
264 | const Int_t bin = (Int_t)x+1;
|
---|
265 | AddBinContent(bin, w);
|
---|
266 | fBinEntries[bin]++;
|
---|
267 | if (fSumw2.fN)
|
---|
268 | fSumw2.fArray[bin] += w*w;
|
---|
269 |
|
---|
270 | if (bin<=0 || bin>fNcells-2)
|
---|
271 | return -1;
|
---|
272 |
|
---|
273 | const Stat_t z = (w > 0 ? w : -w);
|
---|
274 | fTsumw += z;
|
---|
275 | fTsumw2 += z*z;
|
---|
276 | fTsumwx += z*x;
|
---|
277 | fTsumwx2 += z*x*x;
|
---|
278 | return bin;
|
---|
279 | }
|
---|
280 |
|
---|
281 | // ------------------------------------------------------------------------
|
---|
282 | //
|
---|
283 | // Use x and y in millimeters
|
---|
284 | //
|
---|
285 | Int_t MHCamera::Fill(Axis_t x, Axis_t y, Stat_t w)
|
---|
286 | {
|
---|
287 | if (fNcells<=1 || IsFreezed())
|
---|
288 | return -1;
|
---|
289 |
|
---|
290 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
291 | {
|
---|
292 | MHexagon hex((*fGeomCam)[idx]);
|
---|
293 | if (hex.DistanceToPrimitive(x, y)>0)
|
---|
294 | continue;
|
---|
295 |
|
---|
296 | SetUsed(idx);
|
---|
297 | return Fill(idx, w);
|
---|
298 | }
|
---|
299 | return -1;
|
---|
300 | }
|
---|
301 |
|
---|
302 | // ------------------------------------------------------------------------
|
---|
303 | //
|
---|
304 | // Call this if you want to change the display status (displayed or not)
|
---|
305 | // for all pixels. val==0 means that the pixel is not displayed.
|
---|
306 | //
|
---|
307 | void MHCamera::SetUsed(const TArrayC &arr)
|
---|
308 | {
|
---|
309 | if (fNcells-2 != arr.GetSize())
|
---|
310 | {
|
---|
311 | gLog << warn << "WARNING - MHCamera::SetUsed: array size mismatch... ignored." << endl;
|
---|
312 | return;
|
---|
313 | }
|
---|
314 |
|
---|
315 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
316 | arr[idx] ? SetUsed(idx) : ResetUsed(idx);
|
---|
317 | }
|
---|
318 |
|
---|
319 | // ------------------------------------------------------------------------
|
---|
320 | //
|
---|
321 | // Return the mean value of all entries which are used if all=kFALSE and
|
---|
322 | // of all entries if all=kTRUE if sector<0. If sector>=0 only
|
---|
323 | // entries with match the given sector are taken into account.
|
---|
324 | //
|
---|
325 | Stat_t MHCamera::GetMeanSectors(const TArrayI §or, const TArrayI &aidx, Bool_t ball) const
|
---|
326 | {
|
---|
327 | if (fNcells<=1)
|
---|
328 | return 0;
|
---|
329 |
|
---|
330 | Int_t n=0;
|
---|
331 |
|
---|
332 | Stat_t mean = 0;
|
---|
333 | for (int i=0; i<fNcells-2; i++)
|
---|
334 | {
|
---|
335 | if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
|
---|
336 | {
|
---|
337 | if (TestBit(kProfile) && fBinEntries[i+1]==0)
|
---|
338 | continue;
|
---|
339 |
|
---|
340 | mean += TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
|
---|
341 | n++;
|
---|
342 | }
|
---|
343 | }
|
---|
344 |
|
---|
345 | return n==0 ? 0 : mean/n;
|
---|
346 | }
|
---|
347 |
|
---|
348 | // ------------------------------------------------------------------------
|
---|
349 | //
|
---|
350 | // Return the median value of all entries which are used if all=kFALSE and
|
---|
351 | // of all entries if all=kTRUE if sector<0. If sector>=0 only
|
---|
352 | // entries with match the given sector are taken into account.
|
---|
353 | //
|
---|
354 | Stat_t MHCamera::GetMedianSectors(const TArrayI §or, const TArrayI &aidx, Bool_t ball) const
|
---|
355 | {
|
---|
356 | if (fNcells<=1)
|
---|
357 | return 0;
|
---|
358 |
|
---|
359 | TArrayD arr(fNcells-2);
|
---|
360 | Int_t n=0;
|
---|
361 |
|
---|
362 | for (int i=0; i<fNcells-2; i++)
|
---|
363 | {
|
---|
364 | if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
|
---|
365 | {
|
---|
366 | if (TestBit(kProfile) && fBinEntries[i+1]==0)
|
---|
367 | continue;
|
---|
368 |
|
---|
369 | arr[n++] = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
|
---|
370 | }
|
---|
371 | }
|
---|
372 |
|
---|
373 | // return Median of the profile data
|
---|
374 | return TMath::Median(n, arr.GetArray());
|
---|
375 | }
|
---|
376 |
|
---|
377 | // ------------------------------------------------------------------------
|
---|
378 | //
|
---|
379 | // Return the sqrt variance of all entries which are used if all=kFALSE and
|
---|
380 | // of all entries if all=kTRUE if sector<0. If sector>=0 only
|
---|
381 | // entries with match the given sector are taken into account.
|
---|
382 | //
|
---|
383 | Stat_t MHCamera::GetRmsSectors(const TArrayI §or, const TArrayI &aidx, Bool_t ball) const
|
---|
384 | {
|
---|
385 | if (fNcells<=1)
|
---|
386 | return -1;
|
---|
387 |
|
---|
388 | Int_t n=0;
|
---|
389 |
|
---|
390 | Stat_t sum = 0;
|
---|
391 | Stat_t sq = 0;
|
---|
392 | for (int i=0; i<fNcells-2; i++)
|
---|
393 | {
|
---|
394 | if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
|
---|
395 | {
|
---|
396 | if (TestBit(kProfile) && fBinEntries[i+1]==0)
|
---|
397 | continue;
|
---|
398 |
|
---|
399 | const Double_t val = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
|
---|
400 |
|
---|
401 | sum += val;
|
---|
402 | sq += val*val;
|
---|
403 | n++;
|
---|
404 | }
|
---|
405 | }
|
---|
406 |
|
---|
407 | if (n==0)
|
---|
408 | return 0;
|
---|
409 |
|
---|
410 | sum /= n;
|
---|
411 | sq /= n;
|
---|
412 |
|
---|
413 | return TMath::Sqrt(sq-sum*sum);
|
---|
414 | }
|
---|
415 |
|
---|
416 | // ------------------------------------------------------------------------
|
---|
417 | //
|
---|
418 | // Return the median value (divided by MMath::GausProb(1.0)=68.3%) of the
|
---|
419 | // distribution of abs(y[i]-Median). This is my Median equivalent of the RMS.
|
---|
420 | // Return the deviation of all entries which are used if all=kFALSE and
|
---|
421 | // of all entries if all=kTRUE if sector<0. If sector>=0 only
|
---|
422 | // entries with match the given sector are taken into account.
|
---|
423 | //
|
---|
424 | Stat_t MHCamera::GetDevSectors(const TArrayI §or, const TArrayI &aidx, Bool_t ball) const
|
---|
425 | {
|
---|
426 | if (fNcells<=1)
|
---|
427 | return 0;
|
---|
428 |
|
---|
429 | TArrayD arr(fNcells-2);
|
---|
430 | Int_t n=0;
|
---|
431 |
|
---|
432 | for (int i=0; i<fNcells-2; i++)
|
---|
433 | {
|
---|
434 | if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
|
---|
435 | {
|
---|
436 | if (TestBit(kProfile) && fBinEntries[i+1]==0)
|
---|
437 | continue;
|
---|
438 |
|
---|
439 | arr[n++] = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
|
---|
440 | }
|
---|
441 | }
|
---|
442 |
|
---|
443 | // return Median of the profile data
|
---|
444 | return MMath::MedianDev(n, arr.GetArray());
|
---|
445 | }
|
---|
446 |
|
---|
447 | // ------------------------------------------------------------------------
|
---|
448 | //
|
---|
449 | // Return the minimum contents of all pixels (if all is set, otherwise
|
---|
450 | // only of all 'used' pixels), fMinimum if fMinimum set. If sector>=0
|
---|
451 | // only pixels with matching sector number are taken into account.
|
---|
452 | //
|
---|
453 | Double_t MHCamera::GetMinimumSectors(const TArrayI §or, const TArrayI &aidx, Bool_t ball) const
|
---|
454 | {
|
---|
455 | if (fMinimum != -1111)
|
---|
456 | return fMinimum;
|
---|
457 |
|
---|
458 | if (fNcells<=1)
|
---|
459 | return 0;
|
---|
460 |
|
---|
461 | Double_t minimum=FLT_MAX;
|
---|
462 |
|
---|
463 | for (Int_t i=0; i<fNcells-2; i++)
|
---|
464 | {
|
---|
465 | if (TestBit(kProfile) && fBinEntries[i+1]==0)
|
---|
466 | continue;
|
---|
467 |
|
---|
468 | const Double_t val = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
|
---|
469 | if (MatchSector(i, sector, aidx) && (ball || IsUsed(i)) && val<minimum)
|
---|
470 | minimum = val;
|
---|
471 | }
|
---|
472 |
|
---|
473 | return minimum;
|
---|
474 | }
|
---|
475 |
|
---|
476 | // ------------------------------------------------------------------------
|
---|
477 | //
|
---|
478 | // Return the maximum contents of all pixels (if all is set, otherwise
|
---|
479 | // only of all 'used' pixels), fMaximum if fMaximum set. If sector>=0
|
---|
480 | // only pixels with matching sector number are taken into account.
|
---|
481 | //
|
---|
482 | Double_t MHCamera::GetMaximumSectors(const TArrayI §or, const TArrayI &aidx, Bool_t ball) const
|
---|
483 | {
|
---|
484 | if (fMaximum!=-1111)
|
---|
485 | return fMaximum;
|
---|
486 |
|
---|
487 | if (fNcells<=1)
|
---|
488 | return 1;
|
---|
489 |
|
---|
490 | Double_t maximum=-FLT_MAX;
|
---|
491 | for (Int_t i=0; i<fNcells-2; i++)
|
---|
492 | {
|
---|
493 | if (TestBit(kProfile) && fBinEntries[i+1]==0)
|
---|
494 | continue;
|
---|
495 |
|
---|
496 | const Double_t val = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
|
---|
497 | if (MatchSector(i, sector, aidx) && (ball || IsUsed(i)) && val>maximum)
|
---|
498 | maximum = val;
|
---|
499 | }
|
---|
500 |
|
---|
501 | return maximum;
|
---|
502 | }
|
---|
503 |
|
---|
504 | // ------------------------------------------------------------------------
|
---|
505 | //
|
---|
506 | // Call this function to draw the camera layout into your canvas.
|
---|
507 | // Setup a drawing canvas. Add this object and all child objects
|
---|
508 | // (hexagons, etc) to the current pad. If no pad exists a new one is
|
---|
509 | // created. (To access the 'real' pad containing the camera you have
|
---|
510 | // to do a cd(1) in the current layer.
|
---|
511 | //
|
---|
512 | // To draw a camera into its own pad do something like:
|
---|
513 | //
|
---|
514 | // MGeomCamMagic m;
|
---|
515 | // MHCamera *d=new MHCamera(m);
|
---|
516 | //
|
---|
517 | // TCanvas *c = new TCanvas;
|
---|
518 | // c->Divide(2,1);
|
---|
519 | // c->cd(1);
|
---|
520 | //
|
---|
521 | // d->FillRandom();
|
---|
522 | // d->Draw();
|
---|
523 | // d->SetBit(kCanDelete);
|
---|
524 | //
|
---|
525 | // There are several drawing options:
|
---|
526 | // 'hist' Draw as a standard TH1 histogram (value vs. pixel index)
|
---|
527 | // 'box' Draw hexagons which size is in respect to its contents
|
---|
528 | // 'nocol' Leave the 'boxed' hexagons empty
|
---|
529 | // 'pixelindex' Display the pixel index in each pixel
|
---|
530 | // 'sectorindex' Display the sector index in each pixel
|
---|
531 | // 'content' Display the relative content aligned to GetMaximum() and
|
---|
532 | // GeMinimum() ((val-min)/(max-min))
|
---|
533 | // 'proj' Display the y-projection of the histogram
|
---|
534 | // 'pal0' Use Pretty palette
|
---|
535 | // 'pal1' Use Deep Blue Sea palette
|
---|
536 | // 'pal2' Use Inverse Depp Blue Sea palette
|
---|
537 | // 'same' Draw trandparent pixels on top of an existing pad. This
|
---|
538 | // makes it possible to draw the camera image on top of an
|
---|
539 | // existing TH2, but also allows for distorted camera images
|
---|
540 | //
|
---|
541 | void MHCamera::Draw(Option_t *option)
|
---|
542 | {
|
---|
543 | const Bool_t hassame = TString(option).Contains("same", TString::kIgnoreCase) && gPad;
|
---|
544 |
|
---|
545 | // root 3.02:
|
---|
546 | // gPad->SetFixedAspectRatio()
|
---|
547 | const Color_t col = gPad ? gPad->GetFillColor() : 16;
|
---|
548 | TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas("CamDisplay", "Mars Camera Display", 656, 600);
|
---|
549 |
|
---|
550 | if (!hassame)
|
---|
551 | {
|
---|
552 | pad->SetBorderMode(0);
|
---|
553 | pad->SetFillColor(col);
|
---|
554 |
|
---|
555 | //
|
---|
556 | // Create an own pad for the MHCamera-Object which can be
|
---|
557 | // resized in paint to keep the correct aspect ratio
|
---|
558 | //
|
---|
559 | // The margin != 0 is a workaround for a problem in root 4.02/00
|
---|
560 | pad->Divide(1, 1, 1e-10, 1e-10, col);
|
---|
561 | pad->cd(1);
|
---|
562 | gPad->SetBorderMode(0);
|
---|
563 | }
|
---|
564 |
|
---|
565 | AppendPad(option);
|
---|
566 | //fGeomCam->AppendPad();
|
---|
567 |
|
---|
568 | //
|
---|
569 | // Do not change gPad. The user should not see, that Draw
|
---|
570 | // changes gPad...
|
---|
571 | //
|
---|
572 | if (!hassame)
|
---|
573 | pad->cd();
|
---|
574 | }
|
---|
575 |
|
---|
576 | // ------------------------------------------------------------------------
|
---|
577 | //
|
---|
578 | // This is TObject::DrawClone but completely ignores
|
---|
579 | // gROOT->GetSelectedPad(). tbretz had trouble with this in the past.
|
---|
580 | // If this makes trouble please write a bug report.
|
---|
581 | //
|
---|
582 | TObject *MHCamera::DrawClone(Option_t *option) const
|
---|
583 | {
|
---|
584 | // Draw a clone of this object in the current pad
|
---|
585 |
|
---|
586 | //TVirtualPad *pad = gROOT->GetSelectedPad();
|
---|
587 | TVirtualPad *padsav = gPad;
|
---|
588 | //if (pad) pad->cd();
|
---|
589 |
|
---|
590 | TObject *newobj = Clone();
|
---|
591 |
|
---|
592 | if (!newobj)
|
---|
593 | return 0;
|
---|
594 |
|
---|
595 | /*
|
---|
596 | if (pad) {
|
---|
597 | if (strlen(option)) pad->GetListOfPrimitives()->Add(newobj,option);
|
---|
598 | else pad->GetListOfPrimitives()->Add(newobj,GetDrawOption());
|
---|
599 | pad->Modified(kTRUE);
|
---|
600 | pad->Update();
|
---|
601 | if (padsav) padsav->cd();
|
---|
602 | return newobj;
|
---|
603 | }
|
---|
604 | */
|
---|
605 |
|
---|
606 | TString opt(option);
|
---|
607 | opt.ToLower();
|
---|
608 |
|
---|
609 | newobj->Draw(opt.IsNull() ? GetDrawOption() : option);
|
---|
610 |
|
---|
611 | if (padsav)
|
---|
612 | padsav->cd();
|
---|
613 |
|
---|
614 | return newobj;
|
---|
615 | }
|
---|
616 |
|
---|
617 | // ------------------------------------------------------------------------
|
---|
618 | //
|
---|
619 | // Creates a TH1D which contains the projection of the contents of the
|
---|
620 | // MHCamera onto the y-axis. The maximum and minimum are calculated
|
---|
621 | // such that a slighly wider range than (GetMinimum(), GetMaximum()) is
|
---|
622 | // displayed using THLimitsFinder::OptimizeLimits.
|
---|
623 | //
|
---|
624 | // If no name is given the newly allocated histogram is removed from
|
---|
625 | // the current directory calling SetDirectory(0) in any other case
|
---|
626 | // the newly created histogram is removed from the current directory
|
---|
627 | // and added to gROOT such the gROOT->FindObject can find the histogram.
|
---|
628 | //
|
---|
629 | // If the standard name "_py" is given "_py" is appended to the name
|
---|
630 | // of the MHCamera and the corresponding histogram is searched using
|
---|
631 | // gROOT->FindObject and updated with the present projection.
|
---|
632 | //
|
---|
633 | // It is the responsibility of the user to make sure, that the newly
|
---|
634 | // created histogram is freed correctly.
|
---|
635 | //
|
---|
636 | // Currently the new histogram is restrictred to 50 bins.
|
---|
637 | // Maybe a optimal number can be calulated from the number of
|
---|
638 | // bins on the x-axis of the MHCamera?
|
---|
639 | //
|
---|
640 | // The code was taken mainly from TH2::ProjectX such the interface
|
---|
641 | // is more or less the same than to TH2-projections.
|
---|
642 | //
|
---|
643 | // If sector>=0 only entries with matching sector index are taken
|
---|
644 | // into account.
|
---|
645 | //
|
---|
646 | TH1D *MHCamera::ProjectionS(const TArrayI §or, const TArrayI &aidx, const char *name, const Int_t nbins) const
|
---|
647 | {
|
---|
648 |
|
---|
649 | // Create the projection histogram
|
---|
650 | TString pname(name);
|
---|
651 | if (name=="_py")
|
---|
652 | {
|
---|
653 | pname.Prepend(GetName());
|
---|
654 | if (sector.GetSize()>0)
|
---|
655 | {
|
---|
656 | pname += ";";
|
---|
657 | for (int i=0; i<sector.GetSize(); i++)
|
---|
658 | pname += sector[i];
|
---|
659 | }
|
---|
660 | if (aidx.GetSize()>0)
|
---|
661 | {
|
---|
662 | pname += ";";
|
---|
663 | for (int i=0; i<aidx.GetSize(); i++)
|
---|
664 | pname += aidx[i];
|
---|
665 | }
|
---|
666 | }
|
---|
667 |
|
---|
668 | TH1D *h1=0;
|
---|
669 |
|
---|
670 | //check if histogram with identical name exist
|
---|
671 | TObject *h1obj = gROOT->FindObject(pname);
|
---|
672 | if (h1obj && h1obj->InheritsFrom("TH1D")) {
|
---|
673 | h1 = (TH1D*)h1obj;
|
---|
674 | h1->Reset();
|
---|
675 | }
|
---|
676 |
|
---|
677 | if (!h1)
|
---|
678 | {
|
---|
679 | h1 = new TH1D;
|
---|
680 | h1->UseCurrentStyle();
|
---|
681 | h1->SetName(pname);
|
---|
682 | h1->SetTitle(GetTitle());
|
---|
683 | h1->SetDirectory(0);
|
---|
684 | h1->SetXTitle(GetYaxis()->GetTitle());
|
---|
685 | h1->SetYTitle("Counts");
|
---|
686 | //h1->Sumw2();
|
---|
687 | }
|
---|
688 |
|
---|
689 | Double_t min = GetMinimumSectors(sector, aidx);
|
---|
690 | Double_t max = GetMaximumSectors(sector, aidx);
|
---|
691 |
|
---|
692 | if (min==max && max>0)
|
---|
693 | min=0;
|
---|
694 | if (min==max && min<0)
|
---|
695 | max=0;
|
---|
696 |
|
---|
697 | Int_t newbins=0;
|
---|
698 | THLimitsFinder::OptimizeLimits(nbins, newbins, min, max, kFALSE);
|
---|
699 |
|
---|
700 | MBinning bins(nbins, min, max);
|
---|
701 | bins.Apply(*h1);
|
---|
702 |
|
---|
703 | // Fill the projected histogram
|
---|
704 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
705 | if (IsUsed(idx) && MatchSector(idx, sector, aidx))
|
---|
706 | h1->Fill(GetBinContent(idx+1));
|
---|
707 |
|
---|
708 | return h1;
|
---|
709 | }
|
---|
710 |
|
---|
711 | // ------------------------------------------------------------------------
|
---|
712 | //
|
---|
713 | // Creates a TH1D which contains the projection of the contents of the
|
---|
714 | // MHCamera onto the radius from the camera center.
|
---|
715 | // The maximum and minimum are calculated
|
---|
716 | // such that a slighly wider range than (GetMinimum(), GetMaximum()) is
|
---|
717 | // displayed using THLimitsFinder::OptimizeLimits.
|
---|
718 | //
|
---|
719 | // If no name is given the newly allocated histogram is removed from
|
---|
720 | // the current directory calling SetDirectory(0) in any other case
|
---|
721 | // the newly created histogram is removed from the current directory
|
---|
722 | // and added to gROOT such the gROOT->FindObject can find the histogram.
|
---|
723 | //
|
---|
724 | // If the standard name "_rad" is given "_rad" is appended to the name
|
---|
725 | // of the MHCamera and the corresponding histogram is searched using
|
---|
726 | // gROOT->FindObject and updated with the present projection.
|
---|
727 | //
|
---|
728 | // It is the responsibility of the user to make sure, that the newly
|
---|
729 | // created histogram is freed correctly.
|
---|
730 | //
|
---|
731 | // Currently the new histogram is restrictred to 50 bins.
|
---|
732 | // Maybe a optimal number can be calulated from the number of
|
---|
733 | // bins on the x-axis of the MHCamera?
|
---|
734 | //
|
---|
735 | // The code was taken mainly from TH2::ProjectX such the interface
|
---|
736 | // is more or less the same than to TH2-projections.
|
---|
737 | //
|
---|
738 | // If sector>=0 only entries with matching sector index are taken
|
---|
739 | // into account.
|
---|
740 | //
|
---|
741 | TProfile *MHCamera::RadialProfileS(const TArrayI §or, const TArrayI &aidx, const char *name, const Int_t nbins) const
|
---|
742 | {
|
---|
743 | // Create the projection histogram
|
---|
744 | TString pname(name);
|
---|
745 | if (name=="_rad")
|
---|
746 | {
|
---|
747 | pname.Prepend(GetName());
|
---|
748 | if (sector.GetSize()>0)
|
---|
749 | {
|
---|
750 | pname += ";";
|
---|
751 | for (int i=0; i<sector.GetSize(); i++)
|
---|
752 | pname += sector[i];
|
---|
753 | }
|
---|
754 | if (aidx.GetSize()>0)
|
---|
755 | {
|
---|
756 | pname += ";";
|
---|
757 | for (int i=0; i<aidx.GetSize(); i++)
|
---|
758 | pname += aidx[i];
|
---|
759 | }
|
---|
760 | }
|
---|
761 |
|
---|
762 | TProfile *h1=0;
|
---|
763 |
|
---|
764 | //check if histogram with identical name exist
|
---|
765 | TObject *h1obj = gROOT->FindObject(pname);
|
---|
766 | if (h1obj && h1obj->InheritsFrom("TProfile")) {
|
---|
767 | h1 = (TProfile*)h1obj;
|
---|
768 | h1->Reset();
|
---|
769 | }
|
---|
770 |
|
---|
771 | if (!h1)
|
---|
772 | {
|
---|
773 | h1 = new TProfile;
|
---|
774 | h1->UseCurrentStyle();
|
---|
775 | h1->SetName(pname);
|
---|
776 | h1->SetTitle(GetTitle());
|
---|
777 | h1->SetDirectory(0);
|
---|
778 | h1->SetXTitle("Radius from camera center [mm]");
|
---|
779 | h1->SetYTitle(GetYaxis()->GetTitle());
|
---|
780 | }
|
---|
781 |
|
---|
782 | const Double_t m2d = fGeomCam->GetConvMm2Deg();
|
---|
783 |
|
---|
784 | Double_t min = 0.;
|
---|
785 | Double_t max = fGeomCam->GetMaxRadius()*m2d;
|
---|
786 |
|
---|
787 | Int_t newbins=0;
|
---|
788 |
|
---|
789 | THLimitsFinder::OptimizeLimits(nbins, newbins, min, max, kFALSE);
|
---|
790 |
|
---|
791 | MBinning bins(nbins, min, max);
|
---|
792 | bins.Apply(*h1);
|
---|
793 |
|
---|
794 | // Fill the projected histogram
|
---|
795 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
796 | if (IsUsed(idx) && MatchSector(idx, sector, aidx))
|
---|
797 | h1->Fill(TMath::Hypot((*fGeomCam)[idx].GetX(),(*fGeomCam)[idx].GetY())*m2d,
|
---|
798 | GetBinContent(idx+1));
|
---|
799 | return h1;
|
---|
800 | }
|
---|
801 |
|
---|
802 |
|
---|
803 | // ------------------------------------------------------------------------
|
---|
804 | //
|
---|
805 | // Creates a TH1D which contains the projection of the contents of the
|
---|
806 | // MHCamera onto the azimuth angle in the camera.
|
---|
807 | //
|
---|
808 | // If no name is given the newly allocated histogram is removed from
|
---|
809 | // the current directory calling SetDirectory(0) in any other case
|
---|
810 | // the newly created histogram is removed from the current directory
|
---|
811 | // and added to gROOT such the gROOT->FindObject can find the histogram.
|
---|
812 | //
|
---|
813 | // If the standard name "_azi" is given "_azi" is appended to the name
|
---|
814 | // of the MHCamera and the corresponding histogram is searched using
|
---|
815 | // gROOT->FindObject and updated with the present projection.
|
---|
816 | //
|
---|
817 | // It is the responsibility of the user to make sure, that the newly
|
---|
818 | // created histogram is freed correctly.
|
---|
819 | //
|
---|
820 | // Currently the new histogram is restrictred to 60 bins.
|
---|
821 | // Maybe a optimal number can be calulated from the number of
|
---|
822 | // bins on the x-axis of the MHCamera?
|
---|
823 | //
|
---|
824 | // The code was taken mainly from TH2::ProjectX such the interface
|
---|
825 | // is more or less the same than to TH2-projections.
|
---|
826 | //
|
---|
827 | TProfile *MHCamera::AzimuthProfileA(const TArrayI &aidx, const char *name, const Int_t nbins) const
|
---|
828 | {
|
---|
829 | // Create the projection histogram
|
---|
830 | TString pname(name);
|
---|
831 | if (name=="_azi")
|
---|
832 | {
|
---|
833 | pname.Prepend(GetName());
|
---|
834 | if (aidx.GetSize()>0)
|
---|
835 | {
|
---|
836 | pname += ";";
|
---|
837 | for (int i=0; i<aidx.GetSize(); i++)
|
---|
838 | pname += aidx[i];
|
---|
839 | }
|
---|
840 | }
|
---|
841 |
|
---|
842 | TProfile *h1=0;
|
---|
843 |
|
---|
844 | //check if histogram with identical name exist
|
---|
845 | TObject *h1obj = gROOT->FindObject(pname);
|
---|
846 | if (h1obj && h1obj->InheritsFrom("TProfile")) {
|
---|
847 | h1 = (TProfile*)h1obj;
|
---|
848 | h1->Reset();
|
---|
849 | }
|
---|
850 |
|
---|
851 | if (!h1)
|
---|
852 | {
|
---|
853 |
|
---|
854 | h1 = new TProfile;
|
---|
855 | h1->UseCurrentStyle();
|
---|
856 | h1->SetName(pname);
|
---|
857 | h1->SetTitle(GetTitle());
|
---|
858 | h1->SetDirectory(0);
|
---|
859 | h1->SetXTitle("Azimuth in camera [deg]");
|
---|
860 | h1->SetYTitle(GetYaxis()->GetTitle());
|
---|
861 | }
|
---|
862 |
|
---|
863 | //Double_t min = 0;
|
---|
864 | //Double_t max = 360;
|
---|
865 |
|
---|
866 | //Int_t newbins=0;
|
---|
867 | //THLimitsFinder::OptimizeLimits(nbins, newbins, min, max, kFALSE);
|
---|
868 |
|
---|
869 | MBinning bins(nbins, 0, 360);
|
---|
870 | bins.Apply(*h1);
|
---|
871 |
|
---|
872 | // Fill the projected histogram
|
---|
873 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
874 | {
|
---|
875 | if (IsUsed(idx) && MatchSector(idx, TArrayI(), aidx))
|
---|
876 | h1->Fill(TMath::ATan2((*fGeomCam)[idx].GetY(),(*fGeomCam)[idx].GetX())*TMath::RadToDeg()+180,
|
---|
877 | GetPixContent(idx));
|
---|
878 |
|
---|
879 | }
|
---|
880 |
|
---|
881 | return h1;
|
---|
882 | }
|
---|
883 |
|
---|
884 |
|
---|
885 | // ------------------------------------------------------------------------
|
---|
886 | //
|
---|
887 | // Resizes the current pad so that the camera is displayed in its
|
---|
888 | // correct aspect ratio
|
---|
889 | //
|
---|
890 | void MHCamera::SetRange()
|
---|
891 | {
|
---|
892 | const Float_t range = fGeomCam->GetMaxRadius()*1.05;
|
---|
893 |
|
---|
894 | //
|
---|
895 | // Maintain aspect ratio
|
---|
896 | //
|
---|
897 | const float ratio = TestBit(kNoLegend) ? 1 : 1.15;
|
---|
898 |
|
---|
899 | //
|
---|
900 | // Calculate width and height of the current pad in pixels
|
---|
901 | //
|
---|
902 | Float_t w = gPad->GetWw();
|
---|
903 | Float_t h = gPad->GetWh()*ratio;
|
---|
904 |
|
---|
905 | //
|
---|
906 | // This prevents the pad from resizing itself wrongly
|
---|
907 | //
|
---|
908 | if (gPad->GetMother() != gPad)
|
---|
909 | {
|
---|
910 | w *= gPad->GetMother()->GetAbsWNDC();
|
---|
911 | h *= gPad->GetMother()->GetAbsHNDC();
|
---|
912 | }
|
---|
913 |
|
---|
914 | //
|
---|
915 | // Set Range (coordinate system) of pad
|
---|
916 | //
|
---|
917 | gPad->Range(-range, -range, (2*ratio-1)*range, range);
|
---|
918 |
|
---|
919 | //
|
---|
920 | // Resize Pad to given ratio
|
---|
921 | //
|
---|
922 | if (h<w)
|
---|
923 | gPad->SetPad((1.-h/w)/2, 0, (h/w+1.)/2, 1);
|
---|
924 | else
|
---|
925 | gPad->SetPad(0, (1.-w/h)/2, 1, (w/h+1.)/2);
|
---|
926 | }
|
---|
927 |
|
---|
928 | // ------------------------------------------------------------------------
|
---|
929 | //
|
---|
930 | // Updates the pixel colors and paints the pixels
|
---|
931 | //
|
---|
932 | void MHCamera::Update(Bool_t islog, Bool_t isbox, Bool_t iscol, Bool_t issame)
|
---|
933 | {
|
---|
934 | Double_t min = GetMinimum(kFALSE);
|
---|
935 | Double_t max = GetMaximum(kFALSE);
|
---|
936 | if (min==FLT_MAX)
|
---|
937 | {
|
---|
938 | min = 0;
|
---|
939 | max = 1;
|
---|
940 | }
|
---|
941 |
|
---|
942 | if (min==max)
|
---|
943 | max += 1;
|
---|
944 |
|
---|
945 | if (!issame)
|
---|
946 | UpdateLegend(min, max, islog);
|
---|
947 |
|
---|
948 | // Try to estimate the units of the current display. This is only
|
---|
949 | // necessary for 'same' option and allows distorted images of the camera!
|
---|
950 | const Float_t maxr = (1-fGeomCam->GetConvMm2Deg())*fGeomCam->GetMaxRadius()/2;
|
---|
951 | const Float_t conv = !issame ||
|
---|
952 | gPad->GetX1()<-maxr || gPad->GetY1()<-maxr ||
|
---|
953 | gPad->GetX2()> maxr || gPad->GetY2()>maxr ? 1 : fGeomCam->GetConvMm2Deg();
|
---|
954 |
|
---|
955 | MHexagon hex;
|
---|
956 | for (Int_t i=0; i<fNcells-2; i++)
|
---|
957 | {
|
---|
958 | hex.SetFillStyle(issame || (IsTransparent() && !IsUsed(i)) ? 0 : 1001);
|
---|
959 |
|
---|
960 | if (!issame)
|
---|
961 | {
|
---|
962 | const Bool_t isnan = !TMath::Finite(fArray[i+1]);
|
---|
963 | if (!IsUsed(i) || !iscol || isnan)
|
---|
964 | {
|
---|
965 | hex.SetFillColor(10);
|
---|
966 |
|
---|
967 | if (isnan)
|
---|
968 | gLog << warn << "MHCamera::Update: " << GetName() << " <" << GetTitle() << "> - Pixel Index #" << i << " contents is not finite..." << endl;
|
---|
969 | }
|
---|
970 | else
|
---|
971 | hex.SetFillColor(GetColor(GetBinContent(i+1), min, max, islog));
|
---|
972 | }
|
---|
973 |
|
---|
974 | const MGeomPix &pix = (*fGeomCam)[i];
|
---|
975 |
|
---|
976 | Float_t x = pix.GetX()*conv/(fAbberation+1);
|
---|
977 | Float_t y = pix.GetY()*conv/(fAbberation+1);
|
---|
978 | Float_t d = pix.GetD()*conv;
|
---|
979 |
|
---|
980 | if (!isbox)
|
---|
981 | if (IsUsed(i) || !TestBit(kNoUnused))
|
---|
982 | hex.PaintHexagon(x, y, d);
|
---|
983 | else
|
---|
984 | if (IsUsed(i) && TMath::Finite(fArray[i+1]))
|
---|
985 | {
|
---|
986 | Float_t size = d*(GetBinContent(i+1)-min)/(max-min);
|
---|
987 | if (size>d)
|
---|
988 | size=d;
|
---|
989 | hex.PaintHexagon(x, y, size);
|
---|
990 | }
|
---|
991 | }
|
---|
992 | }
|
---|
993 |
|
---|
994 | // ------------------------------------------------------------------------
|
---|
995 | //
|
---|
996 | // Print minimum and maximum
|
---|
997 | //
|
---|
998 | void MHCamera::Print(Option_t *) const
|
---|
999 | {
|
---|
1000 | gLog << all << "Minimum: " << GetMinimum();
|
---|
1001 | if (fMinimum==-1111)
|
---|
1002 | gLog << " <autoscaled>";
|
---|
1003 | gLog << endl;
|
---|
1004 | gLog << "Maximum: " << GetMaximum();
|
---|
1005 | if (fMaximum==-1111)
|
---|
1006 | gLog << " <autoscaled>";
|
---|
1007 | gLog << endl;
|
---|
1008 | }
|
---|
1009 |
|
---|
1010 | // ------------------------------------------------------------------------
|
---|
1011 | //
|
---|
1012 | // Paint the y-axis title
|
---|
1013 | //
|
---|
1014 | void MHCamera::PaintAxisTitle()
|
---|
1015 | {
|
---|
1016 | const Float_t range = fGeomCam->GetMaxRadius()*1.05;
|
---|
1017 | const Float_t w = (1 + 1.5/sqrt((float)(fNcells-2)))*range;
|
---|
1018 |
|
---|
1019 | TLatex *ptitle = new TLatex(w, -.90*range, GetYaxis()->GetTitle());
|
---|
1020 |
|
---|
1021 | ptitle->SetTextSize(0.05);
|
---|
1022 | ptitle->SetTextAlign(21);
|
---|
1023 |
|
---|
1024 | // box with the histogram title
|
---|
1025 | ptitle->SetTextColor(gStyle->GetTitleTextColor());
|
---|
1026 | #if ROOT_VERSION_CODE > ROOT_VERSION(3,05,01)
|
---|
1027 | ptitle->SetTextFont(gStyle->GetTitleFont(""));
|
---|
1028 | #endif
|
---|
1029 | ptitle->Paint();
|
---|
1030 | }
|
---|
1031 |
|
---|
1032 | // ------------------------------------------------------------------------
|
---|
1033 | //
|
---|
1034 | // Paints the camera.
|
---|
1035 | //
|
---|
1036 | void MHCamera::Paint(Option_t *o)
|
---|
1037 | {
|
---|
1038 | if (fNcells<=1)
|
---|
1039 | return;
|
---|
1040 |
|
---|
1041 | TString opt(o);
|
---|
1042 | opt.ToLower();
|
---|
1043 |
|
---|
1044 | if (opt.Contains("hist"))
|
---|
1045 | {
|
---|
1046 | opt.ReplaceAll("hist", "");
|
---|
1047 | opt.ReplaceAll("box", "");
|
---|
1048 | opt.ReplaceAll("pixelindex", "");
|
---|
1049 | opt.ReplaceAll("sectorindex", "");
|
---|
1050 | opt.ReplaceAll("abscontent", "");
|
---|
1051 | opt.ReplaceAll("content", "");
|
---|
1052 | opt.ReplaceAll("proj", "");
|
---|
1053 | opt.ReplaceAll("pal0", "");
|
---|
1054 | opt.ReplaceAll("pal1", "");
|
---|
1055 | opt.ReplaceAll("pal2", "");
|
---|
1056 | opt.ReplaceAll("nopal", "");
|
---|
1057 | TH1D::Paint(opt);
|
---|
1058 | return;
|
---|
1059 | }
|
---|
1060 |
|
---|
1061 | if (opt.Contains("proj"))
|
---|
1062 | {
|
---|
1063 | opt.ReplaceAll("proj", "");
|
---|
1064 | Projection(GetName())->Paint(opt);
|
---|
1065 | return;
|
---|
1066 | }
|
---|
1067 |
|
---|
1068 | const Bool_t hassame = opt.Contains("same");
|
---|
1069 | const Bool_t hasbox = opt.Contains("box");
|
---|
1070 | const Bool_t hascol = hasbox ? !opt.Contains("nocol") : kTRUE;
|
---|
1071 |
|
---|
1072 | if (!hassame)
|
---|
1073 | {
|
---|
1074 | gPad->Clear();
|
---|
1075 |
|
---|
1076 | // Maintain aspect ratio
|
---|
1077 | SetRange();
|
---|
1078 |
|
---|
1079 | if (GetPainter())
|
---|
1080 | {
|
---|
1081 | // Paint statistics
|
---|
1082 | if (!TestBit(TH1::kNoStats))
|
---|
1083 | fPainter->PaintStat(gStyle->GetOptStat(), NULL);
|
---|
1084 |
|
---|
1085 | // Paint primitives (pixels, color legend, photons, ...)
|
---|
1086 | if (fPainter->InheritsFrom(THistPainter::Class()))
|
---|
1087 | {
|
---|
1088 | static_cast<THistPainter*>(fPainter)->MakeChopt("");
|
---|
1089 | static_cast<THistPainter*>(fPainter)->PaintTitle();
|
---|
1090 | }
|
---|
1091 | }
|
---|
1092 | }
|
---|
1093 |
|
---|
1094 | const Bool_t pal1 = opt.Contains("pal1");
|
---|
1095 | const Bool_t pal2 = opt.Contains("pal2");
|
---|
1096 | const Bool_t nopal = opt.Contains("nopal");
|
---|
1097 |
|
---|
1098 | if (!pal1 && !pal2 && !nopal)
|
---|
1099 | SetPrettyPalette();
|
---|
1100 |
|
---|
1101 | if (pal1)
|
---|
1102 | SetDeepBlueSeaPalette();
|
---|
1103 |
|
---|
1104 | if (pal2)
|
---|
1105 | SetInvDeepBlueSeaPalette();
|
---|
1106 |
|
---|
1107 | // Update Contents of the pixels and paint legend
|
---|
1108 | Update(gPad->GetLogy(), hasbox, hascol, hassame);
|
---|
1109 |
|
---|
1110 | if (!hassame)
|
---|
1111 | PaintAxisTitle();
|
---|
1112 |
|
---|
1113 | if (opt.Contains("pixelindex"))
|
---|
1114 | {
|
---|
1115 | PaintIndices(0);
|
---|
1116 | return;
|
---|
1117 | }
|
---|
1118 | if (opt.Contains("sectorindex"))
|
---|
1119 | {
|
---|
1120 | PaintIndices(1);
|
---|
1121 | return;
|
---|
1122 | }
|
---|
1123 | if (opt.Contains("abscontent"))
|
---|
1124 | {
|
---|
1125 | PaintIndices(3);
|
---|
1126 | return;
|
---|
1127 | }
|
---|
1128 | if (opt.Contains("content"))
|
---|
1129 | {
|
---|
1130 | PaintIndices(2);
|
---|
1131 | return;
|
---|
1132 | }
|
---|
1133 | if (opt.Contains("pixelentries"))
|
---|
1134 | {
|
---|
1135 | PaintIndices(4);
|
---|
1136 | return;
|
---|
1137 | }
|
---|
1138 | }
|
---|
1139 |
|
---|
1140 | void MHCamera::SetDrawOption(Option_t *option)
|
---|
1141 | {
|
---|
1142 | // This is a workaround. For some reason MHCamera is
|
---|
1143 | // stored in a TObjLink instead of a TObjOptLink
|
---|
1144 | if (!option || !gPad)
|
---|
1145 | return;
|
---|
1146 |
|
---|
1147 | TListIter next(gPad->GetListOfPrimitives());
|
---|
1148 | delete gPad->FindObject("Tframe");
|
---|
1149 | TObject *obj;
|
---|
1150 | while ((obj = next()))
|
---|
1151 | if (obj == this && (TString)next.GetOption()!=(TString)option)
|
---|
1152 | {
|
---|
1153 | gPad->GetListOfPrimitives()->Remove(this);
|
---|
1154 | gPad->GetListOfPrimitives()->AddFirst(this, option);
|
---|
1155 | return;
|
---|
1156 | }
|
---|
1157 | }
|
---|
1158 |
|
---|
1159 | // ------------------------------------------------------------------------
|
---|
1160 | //
|
---|
1161 | // With this function you can change the color palette. For more
|
---|
1162 | // information see TStyle::SetPalette. Only palettes with 50 colors
|
---|
1163 | // are allowed.
|
---|
1164 | // In addition you can use SetPalette(52, 0) to create an inverse
|
---|
1165 | // deep blue sea palette
|
---|
1166 | //
|
---|
1167 | void MHCamera::SetPalette(Int_t ncolors, Int_t *colors)
|
---|
1168 | {
|
---|
1169 | //
|
---|
1170 | // If not enough colors are specified skip this.
|
---|
1171 | //
|
---|
1172 | if (ncolors>1 && ncolors<50)
|
---|
1173 | {
|
---|
1174 | gLog << err << "MHCamera::SetPalette: Only default palettes with 50 colors are allowed... ignored." << endl;
|
---|
1175 | return;
|
---|
1176 | }
|
---|
1177 |
|
---|
1178 | //
|
---|
1179 | // If ncolors==52 create a reversed deep blue sea palette
|
---|
1180 | //
|
---|
1181 | if (ncolors==52)
|
---|
1182 | {
|
---|
1183 | gStyle->SetPalette(51, NULL);
|
---|
1184 | TArrayI c(kItemsLegend);
|
---|
1185 | for (int i=0; i<kItemsLegend; i++)
|
---|
1186 | c[kItemsLegend-i-1] = gStyle->GetColorPalette(i);
|
---|
1187 | gStyle->SetPalette(kItemsLegend, c.GetArray());
|
---|
1188 | }
|
---|
1189 | else
|
---|
1190 | gStyle->SetPalette(ncolors, colors);
|
---|
1191 | }
|
---|
1192 |
|
---|
1193 |
|
---|
1194 | // ------------------------------------------------------------------------
|
---|
1195 | //
|
---|
1196 | // Changes the palette of the displayed camera histogram.
|
---|
1197 | //
|
---|
1198 | // Change to the right pad first - otherwise GetDrawOption() might fail.
|
---|
1199 | //
|
---|
1200 | void MHCamera::SetPrettyPalette()
|
---|
1201 | {
|
---|
1202 | TString opt(GetDrawOption());
|
---|
1203 |
|
---|
1204 | if (!opt.Contains("hist", TString::kIgnoreCase))
|
---|
1205 | SetPalette(1, 0);
|
---|
1206 |
|
---|
1207 | opt.ReplaceAll("pal1", "");
|
---|
1208 | opt.ReplaceAll("pal2", "");
|
---|
1209 |
|
---|
1210 | SetDrawOption(opt);
|
---|
1211 | }
|
---|
1212 |
|
---|
1213 | // ------------------------------------------------------------------------
|
---|
1214 | //
|
---|
1215 | // Changes the palette of the displayed camera histogram.
|
---|
1216 | //
|
---|
1217 | // Change to the right pad first - otherwise GetDrawOption() might fail.
|
---|
1218 | //
|
---|
1219 | void MHCamera::SetDeepBlueSeaPalette()
|
---|
1220 | {
|
---|
1221 | TString opt(GetDrawOption());
|
---|
1222 |
|
---|
1223 | if (!opt.Contains("hist", TString::kIgnoreCase))
|
---|
1224 | SetPalette(51, 0);
|
---|
1225 |
|
---|
1226 | opt.ReplaceAll("pal1", "");
|
---|
1227 | opt.ReplaceAll("pal2", "");
|
---|
1228 | opt += "pal1";
|
---|
1229 |
|
---|
1230 | SetDrawOption(opt);
|
---|
1231 | }
|
---|
1232 |
|
---|
1233 | // ------------------------------------------------------------------------
|
---|
1234 | //
|
---|
1235 | // Changes the palette of the displayed camera histogram.
|
---|
1236 | //
|
---|
1237 | // Change to the right pad first - otherwise GetDrawOption() might fail.
|
---|
1238 | //
|
---|
1239 | void MHCamera::SetInvDeepBlueSeaPalette()
|
---|
1240 | {
|
---|
1241 | TString opt(GetDrawOption());
|
---|
1242 |
|
---|
1243 | if (!opt.Contains("hist", TString::kIgnoreCase))
|
---|
1244 | SetPalette(52, 0);
|
---|
1245 |
|
---|
1246 | opt.ReplaceAll("pal1", "");
|
---|
1247 | opt.ReplaceAll("pal2", "");
|
---|
1248 | opt += "pal2";
|
---|
1249 |
|
---|
1250 | SetDrawOption(opt);
|
---|
1251 | }
|
---|
1252 |
|
---|
1253 | // ------------------------------------------------------------------------
|
---|
1254 | //
|
---|
1255 | // Paint indices (as text) inside the pixels. Depending of the type-
|
---|
1256 | // argument we paint:
|
---|
1257 | // 0: pixel number
|
---|
1258 | // 1: sector number
|
---|
1259 | // 2: content
|
---|
1260 | //
|
---|
1261 | void MHCamera::PaintIndices(Int_t type)
|
---|
1262 | {
|
---|
1263 | if (fNcells<=1)
|
---|
1264 | return;
|
---|
1265 |
|
---|
1266 | const Double_t min = GetMinimum();
|
---|
1267 | const Double_t max = GetMaximum();
|
---|
1268 |
|
---|
1269 | if (type==2 && max==min)
|
---|
1270 | return;
|
---|
1271 |
|
---|
1272 | TText txt;
|
---|
1273 | txt.SetTextFont(122);
|
---|
1274 | txt.SetTextAlign(22); // centered/centered
|
---|
1275 |
|
---|
1276 | for (Int_t i=0; i<fNcells-2; i++)
|
---|
1277 | {
|
---|
1278 | const MGeomPix &h = (*fGeomCam)[i];
|
---|
1279 |
|
---|
1280 | TString num;
|
---|
1281 | switch (type)
|
---|
1282 | {
|
---|
1283 | case 0: num += i; break;
|
---|
1284 | case 1: num += h.GetSector(); break;
|
---|
1285 | case 2: num += TMath::Nint((fArray[i+1]-min)/(max-min)); break;
|
---|
1286 | case 3: num += TMath::Nint(fArray[i+1]); break;
|
---|
1287 | case 4: num += fBinEntries[i+1]; break;
|
---|
1288 | }
|
---|
1289 |
|
---|
1290 | // FIXME: Should depend on the color of the pixel...
|
---|
1291 | //(GetColor(GetBinContent(i+1), min, max, 0));
|
---|
1292 | txt.SetTextColor(kRed);
|
---|
1293 | txt.SetTextSize(0.3*h.GetD()/fGeomCam->GetMaxRadius()/1.05);
|
---|
1294 | txt.PaintText(h.GetX(), h.GetY(), num);
|
---|
1295 | }
|
---|
1296 | }
|
---|
1297 |
|
---|
1298 | // ------------------------------------------------------------------------
|
---|
1299 | //
|
---|
1300 | // Call this function to add a MCamEvent on top of the present contents.
|
---|
1301 | //
|
---|
1302 | void MHCamera::AddCamContent(const MCamEvent &event, Int_t type)
|
---|
1303 | {
|
---|
1304 | if (fNcells<=1 || IsFreezed())
|
---|
1305 | return;
|
---|
1306 |
|
---|
1307 | // FIXME: Security check missing!
|
---|
1308 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1309 | {
|
---|
1310 | Double_t val=0;
|
---|
1311 | if (event.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
|
---|
1312 | {
|
---|
1313 | SetUsed(idx);
|
---|
1314 | Fill(idx, val); // FIXME: Slow!
|
---|
1315 | }
|
---|
1316 | }
|
---|
1317 | fEntries++;
|
---|
1318 | }
|
---|
1319 |
|
---|
1320 | // ------------------------------------------------------------------------
|
---|
1321 | //
|
---|
1322 | // Call this function to add a MCamEvent on top of the present contents.
|
---|
1323 | //
|
---|
1324 | void MHCamera::SetCamError(const MCamEvent &evt, Int_t type)
|
---|
1325 | {
|
---|
1326 |
|
---|
1327 | if (fNcells<=1 || IsFreezed())
|
---|
1328 | return;
|
---|
1329 |
|
---|
1330 | // FIXME: Security check missing!
|
---|
1331 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1332 | {
|
---|
1333 | Double_t val=0;
|
---|
1334 | if (evt.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
|
---|
1335 | SetUsed(idx);
|
---|
1336 |
|
---|
1337 | SetBinError(idx+1, val); // FIXME: Slow!
|
---|
1338 | }
|
---|
1339 | }
|
---|
1340 |
|
---|
1341 | Stat_t MHCamera::GetBinContent(Int_t bin) const
|
---|
1342 | {
|
---|
1343 | if (fBuffer) ((TH1D*)this)->BufferEmpty();
|
---|
1344 | if (bin < 0) bin = 0;
|
---|
1345 | if (bin >= fNcells) bin = fNcells-1;
|
---|
1346 | if (!fArray) return 0;
|
---|
1347 |
|
---|
1348 | if (!TestBit(kProfile))
|
---|
1349 | return Stat_t (fArray[bin]);
|
---|
1350 |
|
---|
1351 | if (fBinEntries.fArray[bin] == 0) return 0;
|
---|
1352 | return fArray[bin]/fBinEntries.fArray[bin];
|
---|
1353 | }
|
---|
1354 |
|
---|
1355 | // ------------------------------------------------------------------------
|
---|
1356 | //
|
---|
1357 | // In the case the kProfile flag is set the spread of the bin is returned.
|
---|
1358 | // If you want to have the mean error instead set the kErrorMean bit via
|
---|
1359 | // SetBit(kErrorMean) first.
|
---|
1360 | //
|
---|
1361 | Stat_t MHCamera::GetBinError(Int_t bin) const
|
---|
1362 | {
|
---|
1363 | if (!TestBit(kProfile))
|
---|
1364 | return TH1D::GetBinError(bin);
|
---|
1365 |
|
---|
1366 | const UInt_t n = (UInt_t)fBinEntries[bin];
|
---|
1367 |
|
---|
1368 | if (n==0)
|
---|
1369 | return 0;
|
---|
1370 |
|
---|
1371 | const Double_t sqr = fSumw2.fArray[bin] / n;
|
---|
1372 | const Double_t val = fArray[bin] / n;
|
---|
1373 |
|
---|
1374 | const Double_t spread = sqr>val*val ? TMath::Sqrt(sqr - val*val) : 0;
|
---|
1375 |
|
---|
1376 | return TestBit(kErrorMean) ? spread/TMath::Sqrt(n) : spread;
|
---|
1377 |
|
---|
1378 | /*
|
---|
1379 | Double_t rc = 0;
|
---|
1380 | if (TestBit(kSqrtVariance) && GetEntries()>0) // error on the mean
|
---|
1381 | {
|
---|
1382 | const Double_t error = fSumw2.fArray[bin]/GetEntries();
|
---|
1383 | const Double_t val = fArray[bin]/GetEntries();
|
---|
1384 | rc = val*val>error ? 0 : TMath::Sqrt(error - val*val);
|
---|
1385 | }
|
---|
1386 | else
|
---|
1387 | rc = TH1D::GetBinError(bin);
|
---|
1388 |
|
---|
1389 | return Profile(rc);*/
|
---|
1390 | }
|
---|
1391 |
|
---|
1392 | // ------------------------------------------------------------------------
|
---|
1393 | //
|
---|
1394 | // Call this function to add a MHCamera on top of the present contents.
|
---|
1395 | // Type:
|
---|
1396 | // 0) bin content
|
---|
1397 | // 1) errors
|
---|
1398 | // 2) rel. errors
|
---|
1399 | //
|
---|
1400 | void MHCamera::AddCamContent(const MHCamera &d, Int_t type)
|
---|
1401 | {
|
---|
1402 | if (fNcells!=d.fNcells || IsFreezed())
|
---|
1403 | return;
|
---|
1404 |
|
---|
1405 | // FIXME: Security check missing!
|
---|
1406 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1407 | if (d.IsUsed(idx))
|
---|
1408 | SetUsed(idx);
|
---|
1409 |
|
---|
1410 | switch (type)
|
---|
1411 | {
|
---|
1412 | case 1:
|
---|
1413 | // Under-/Overflow bins not handled!
|
---|
1414 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1415 | if (d.IsUsed(idx))
|
---|
1416 | Fill(idx, d.GetBinError(idx+1));
|
---|
1417 | fEntries++;
|
---|
1418 | break;
|
---|
1419 | case 2:
|
---|
1420 | // Under-/Overflow bins not handled!
|
---|
1421 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1422 | if (d.GetBinContent(idx+1)!=0 && d.IsUsed(idx))
|
---|
1423 | Fill(idx, TMath::Abs(d.GetBinError(idx+1)/d.GetBinContent(idx+1)));
|
---|
1424 | fEntries++;
|
---|
1425 | break;
|
---|
1426 | default:
|
---|
1427 | if (TestBit(kProfile)!=d.TestBit(kProfile))
|
---|
1428 | gLog << warn << "WARNING - You have tried to call AddCamContent for two different kind of histograms (kProfile set or not)." << endl;
|
---|
1429 |
|
---|
1430 | // environment
|
---|
1431 | fEntries += d.fEntries;
|
---|
1432 | fTsumw += d.fTsumw;
|
---|
1433 | fTsumw2 += d.fTsumw2;
|
---|
1434 | fTsumwx += d.fTsumwx;
|
---|
1435 | fTsumwx2 += d.fTsumwx2;
|
---|
1436 | // Bin contents
|
---|
1437 | for (Int_t idx=1; idx<fNcells-1; idx++)
|
---|
1438 | {
|
---|
1439 | if (!d.IsUsed(idx-1))
|
---|
1440 | continue;
|
---|
1441 |
|
---|
1442 | fArray[idx] += d.fArray[idx];
|
---|
1443 | fBinEntries[idx] += d.fBinEntries[idx];
|
---|
1444 | fSumw2.fArray[idx] += d.fSumw2.fArray[idx];
|
---|
1445 | }
|
---|
1446 | // Underflow bin
|
---|
1447 | fArray[0] += d.fArray[0];
|
---|
1448 | fBinEntries[0] += d.fBinEntries[0];
|
---|
1449 | fSumw2.fArray[0] += d.fSumw2.fArray[0];
|
---|
1450 | // Overflow bin
|
---|
1451 | fArray[fNcells-1] += d.fArray[fNcells-1];
|
---|
1452 | fBinEntries[fNcells-1] += d.fBinEntries[fNcells-1];
|
---|
1453 | fSumw2.fArray[fNcells-1] += d.fSumw2.fArray[fNcells-1];
|
---|
1454 | break;
|
---|
1455 | /* default:
|
---|
1456 | if (TestBit(kProfile)!=d.TestBit(kProfile))
|
---|
1457 | gLog << warn << "WARNING - You have tried to call AddCamContent for two different kind of histograms (kProfile set or not)." << endl;
|
---|
1458 |
|
---|
1459 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1460 | Fill(idx, d.GetBinContent(idx+1));
|
---|
1461 | break;*/
|
---|
1462 | }
|
---|
1463 | fEntries++;
|
---|
1464 | }
|
---|
1465 |
|
---|
1466 | // ------------------------------------------------------------------------
|
---|
1467 | //
|
---|
1468 | // Call this function to add a TArrayD on top of the present contents.
|
---|
1469 | //
|
---|
1470 | void MHCamera::AddCamContent(const TArrayD &event, const TArrayC *used)
|
---|
1471 | {
|
---|
1472 | if (event.GetSize()!=fNcells-2 || IsFreezed())
|
---|
1473 | return;
|
---|
1474 |
|
---|
1475 | if (used && used->GetSize()!=fNcells-2)
|
---|
1476 | return;
|
---|
1477 |
|
---|
1478 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1479 | {
|
---|
1480 | Fill(idx, event[idx]); // FIXME: Slow!
|
---|
1481 |
|
---|
1482 | if (!used || (*used)[idx])
|
---|
1483 | SetUsed(idx);
|
---|
1484 | }
|
---|
1485 | fEntries++;
|
---|
1486 | }
|
---|
1487 |
|
---|
1488 | // ------------------------------------------------------------------------
|
---|
1489 | //
|
---|
1490 | // Call this function to add a MArrayD on top of the present contents.
|
---|
1491 | //
|
---|
1492 | void MHCamera::AddCamContent(const MArrayD &event, const TArrayC *used)
|
---|
1493 | {
|
---|
1494 | if (event.GetSize()!=(UInt_t)(fNcells-2) || IsFreezed())
|
---|
1495 | return;
|
---|
1496 |
|
---|
1497 | if (used && used->GetSize()!=fNcells-2)
|
---|
1498 | return;
|
---|
1499 |
|
---|
1500 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1501 | {
|
---|
1502 | Fill(idx, event[idx]); // FIXME: Slow!
|
---|
1503 |
|
---|
1504 | if (!used || (*used)[idx])
|
---|
1505 | SetUsed(idx);
|
---|
1506 | }
|
---|
1507 | fEntries++;
|
---|
1508 | }
|
---|
1509 |
|
---|
1510 | // ------------------------------------------------------------------------
|
---|
1511 | //
|
---|
1512 | // Call this function to add a MCamEvent on top of the present contents.
|
---|
1513 | // 1 is added to each pixel if the contents of MCamEvent>threshold (in case isabove is set to kTRUE == default)
|
---|
1514 | // 1 is added to each pixel if the contents of MCamEvent<threshold (in case isabove is set to kFALSE)
|
---|
1515 | //
|
---|
1516 | // in unused pixel is not counted if it didn't fullfill the condition.
|
---|
1517 | //
|
---|
1518 | void MHCamera::CntCamContent(const MCamEvent &event, Double_t threshold, Int_t type, Bool_t isabove)
|
---|
1519 | {
|
---|
1520 | if (fNcells<=1 || IsFreezed())
|
---|
1521 | return;
|
---|
1522 |
|
---|
1523 | // FIXME: Security check missing!
|
---|
1524 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1525 | {
|
---|
1526 | Double_t val=threshold;
|
---|
1527 | const Bool_t rc = event.GetPixelContent(val, idx, *fGeomCam, type);
|
---|
1528 | if (rc)
|
---|
1529 | SetUsed(idx);
|
---|
1530 |
|
---|
1531 | const Bool_t cond =
|
---|
1532 | ( isabove && val>threshold) ||
|
---|
1533 | (!isabove && val<threshold);
|
---|
1534 |
|
---|
1535 | Fill(idx, rc && cond ? 1 : 0);
|
---|
1536 | }
|
---|
1537 | fEntries++;
|
---|
1538 | }
|
---|
1539 |
|
---|
1540 | // ------------------------------------------------------------------------
|
---|
1541 | //
|
---|
1542 | // Call this function to add a MCamEvent on top of the present contents.
|
---|
1543 | // - the contents of the pixels in event are added to each pixel
|
---|
1544 | // if the pixel of thresevt<threshold (in case isabove is set
|
---|
1545 | // to kTRUE == default)
|
---|
1546 | // - the contents of the pixels in event are added to each pixel
|
---|
1547 | // if the pixel of thresevt<threshold (in case isabove is set
|
---|
1548 | // to kFALSE)
|
---|
1549 | //
|
---|
1550 | // in unused pixel is not counted if it didn't fullfill the condition.
|
---|
1551 | //
|
---|
1552 | void MHCamera::CntCamContent(const MCamEvent &event, Int_t type1, const MCamEvent &thresevt, Int_t type2, Double_t threshold, Bool_t isabove)
|
---|
1553 | {
|
---|
1554 | if (fNcells<=1 || IsFreezed())
|
---|
1555 | return;
|
---|
1556 |
|
---|
1557 | // FIXME: Security check missing!
|
---|
1558 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1559 | {
|
---|
1560 | Double_t th=0;
|
---|
1561 | if (!thresevt.GetPixelContent(th, idx, *fGeomCam, type2))
|
---|
1562 | continue;
|
---|
1563 |
|
---|
1564 | if ((isabove && th>threshold) || (!isabove && th<threshold))
|
---|
1565 | continue;
|
---|
1566 |
|
---|
1567 | Double_t val=th;
|
---|
1568 | if (event.GetPixelContent(val, idx, *fGeomCam, type1))
|
---|
1569 | {
|
---|
1570 | SetUsed(idx);
|
---|
1571 | Fill(idx, val);
|
---|
1572 | }
|
---|
1573 | }
|
---|
1574 | fEntries++;
|
---|
1575 | }
|
---|
1576 |
|
---|
1577 | // ------------------------------------------------------------------------
|
---|
1578 | //
|
---|
1579 | // Call this function to add a MCamEvent on top of the present contents.
|
---|
1580 | // 1 is added to each pixel if the contents of MCamEvent>threshold (in case isabove is set to kTRUE == default)
|
---|
1581 | // 1 is added to each pixel if the contents of MCamEvent<threshold (in case isabove is set to kFALSE)
|
---|
1582 | //
|
---|
1583 | // in unused pixel is not counted if it didn't fullfill the condition.
|
---|
1584 | //
|
---|
1585 | void MHCamera::CntCamContent(const MCamEvent &event, TArrayD threshold, Int_t type, Bool_t isabove)
|
---|
1586 | {
|
---|
1587 | if (fNcells<=1 || IsFreezed())
|
---|
1588 | return;
|
---|
1589 |
|
---|
1590 | // FIXME: Security check missing!
|
---|
1591 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1592 | {
|
---|
1593 | Double_t val=threshold[idx];
|
---|
1594 | if (event.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
|
---|
1595 | {
|
---|
1596 | SetUsed(idx);
|
---|
1597 |
|
---|
1598 | if (val>threshold[idx] && isabove)
|
---|
1599 | Fill(idx);
|
---|
1600 | if (val<threshold[idx] && !isabove)
|
---|
1601 | Fill(idx);
|
---|
1602 | }
|
---|
1603 | }
|
---|
1604 | fEntries++;
|
---|
1605 | }
|
---|
1606 |
|
---|
1607 | // ------------------------------------------------------------------------
|
---|
1608 | //
|
---|
1609 | // Call this function to add a TArrayD on top of the present contents.
|
---|
1610 | // 1 is added to each pixel if the contents of MCamEvent>threshold
|
---|
1611 | //
|
---|
1612 | void MHCamera::CntCamContent(const TArrayD &event, Double_t threshold, Bool_t ispos)
|
---|
1613 | {
|
---|
1614 | if (event.GetSize()!=fNcells-2 || IsFreezed())
|
---|
1615 | return;
|
---|
1616 |
|
---|
1617 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1618 | {
|
---|
1619 | if (event[idx]>threshold)
|
---|
1620 | Fill(idx);
|
---|
1621 |
|
---|
1622 | if (!ispos || fArray[idx+1]>0)
|
---|
1623 | SetUsed(idx);
|
---|
1624 | }
|
---|
1625 | fEntries++;
|
---|
1626 | }
|
---|
1627 |
|
---|
1628 | // ------------------------------------------------------------------------
|
---|
1629 | //
|
---|
1630 | // Fill the pixels with random contents.
|
---|
1631 | //
|
---|
1632 | void MHCamera::FillRandom()
|
---|
1633 | {
|
---|
1634 | if (fNcells<=1 || IsFreezed())
|
---|
1635 | return;
|
---|
1636 |
|
---|
1637 | Reset();
|
---|
1638 |
|
---|
1639 | // FIXME: Security check missing!
|
---|
1640 | for (Int_t idx=0; idx<fNcells-2; idx++)
|
---|
1641 | {
|
---|
1642 | Fill(idx, gRandom->Uniform()*fGeomCam->GetPixRatio(idx));
|
---|
1643 | SetUsed(idx);
|
---|
1644 | }
|
---|
1645 | fEntries=1;
|
---|
1646 | }
|
---|
1647 |
|
---|
1648 |
|
---|
1649 | // ------------------------------------------------------------------------
|
---|
1650 | //
|
---|
1651 | // The array must be in increasing order, eg: 2.5, 3.7, 4.9
|
---|
1652 | // The values in each bin are replaced by the interval in which the value
|
---|
1653 | // fits. In the example we have four intervals
|
---|
1654 | // (<2.5, 2.5-3.7, 3.7-4.9, >4.9). Maximum and minimum are set
|
---|
1655 | // accordingly.
|
---|
1656 | //
|
---|
1657 | void MHCamera::SetLevels(const TArrayF &arr)
|
---|
1658 | {
|
---|
1659 | if (fNcells<=1)
|
---|
1660 | return;
|
---|
1661 |
|
---|
1662 | for (Int_t i=0; i<fNcells-2; i++)
|
---|
1663 | {
|
---|
1664 | if (!IsUsed(i))
|
---|
1665 | continue;
|
---|
1666 |
|
---|
1667 | Int_t j = arr.GetSize();
|
---|
1668 | while (j && fArray[i+1]<arr[j-1])
|
---|
1669 | j--;
|
---|
1670 |
|
---|
1671 | fArray[i+1] = j;
|
---|
1672 | }
|
---|
1673 | SetMaximum(arr.GetSize());
|
---|
1674 | SetMinimum(0);
|
---|
1675 | }
|
---|
1676 |
|
---|
1677 | // ------------------------------------------------------------------------
|
---|
1678 | //
|
---|
1679 | // Reset the all pixel colors to a default value
|
---|
1680 | //
|
---|
1681 | void MHCamera::Reset(Option_t *opt)
|
---|
1682 | {
|
---|
1683 | if (fNcells<=1 || IsFreezed())
|
---|
1684 | return;
|
---|
1685 |
|
---|
1686 | TH1::Reset(opt);
|
---|
1687 |
|
---|
1688 | fUsed.Reset();
|
---|
1689 | fBinEntries.Reset();
|
---|
1690 |
|
---|
1691 | for (Int_t i=0; i<fNcells; i++)
|
---|
1692 | fArray[i] = 0;
|
---|
1693 | }
|
---|
1694 |
|
---|
1695 | // ------------------------------------------------------------------------
|
---|
1696 | //
|
---|
1697 | // Here we calculate the color index for the current value.
|
---|
1698 | // The color index is defined with the class TStyle and the
|
---|
1699 | // Color palette inside. We use the command gStyle->SetPalette(1,0)
|
---|
1700 | // for the display. So we have to convert the value "wert" into
|
---|
1701 | // a color index that fits the color palette.
|
---|
1702 | // The range of the color palette is defined by the values fMinPhe
|
---|
1703 | // and fMaxRange. Between this values we have 50 color index, starting
|
---|
1704 | // with 0 up to 49.
|
---|
1705 | //
|
---|
1706 | Int_t MHCamera::GetColor(Float_t val, Float_t min, Float_t max, Bool_t islog)
|
---|
1707 | {
|
---|
1708 | //
|
---|
1709 | // first treat the over- and under-flows
|
---|
1710 | //
|
---|
1711 | const Int_t maxcolidx = kItemsLegend-1;
|
---|
1712 |
|
---|
1713 | if (!TMath::Finite(val)) // FIXME: gLog!
|
---|
1714 | return maxcolidx/2;
|
---|
1715 |
|
---|
1716 | if (val >= max)
|
---|
1717 | return gStyle->GetColorPalette(maxcolidx);
|
---|
1718 |
|
---|
1719 | if (val <= min)
|
---|
1720 | return gStyle->GetColorPalette(0);
|
---|
1721 |
|
---|
1722 | //
|
---|
1723 | // calculate the color index
|
---|
1724 | //
|
---|
1725 | Float_t ratio;
|
---|
1726 | if (islog && min>0)
|
---|
1727 | ratio = log10(val/min) / log10(max/min);
|
---|
1728 | else
|
---|
1729 | ratio = (val-min) / (max-min);
|
---|
1730 |
|
---|
1731 | const Int_t colidx = (Int_t)(ratio*maxcolidx + .5);
|
---|
1732 | return gStyle->GetColorPalette(colidx);
|
---|
1733 | }
|
---|
1734 |
|
---|
1735 | TPaveStats *MHCamera::GetStatisticBox()
|
---|
1736 | {
|
---|
1737 | TObject *obj = 0;
|
---|
1738 |
|
---|
1739 | TIter Next(fFunctions);
|
---|
1740 | while ((obj = Next()))
|
---|
1741 | if (obj->InheritsFrom(TPaveStats::Class()))
|
---|
1742 | return static_cast<TPaveStats*>(obj);
|
---|
1743 |
|
---|
1744 | return NULL;
|
---|
1745 | }
|
---|
1746 |
|
---|
1747 | // ------------------------------------------------------------------------
|
---|
1748 | //
|
---|
1749 | // Change the text on the legend according to the range of the Display
|
---|
1750 | //
|
---|
1751 | void MHCamera::UpdateLegend(Float_t min, Float_t max, Bool_t islog)
|
---|
1752 | {
|
---|
1753 | const Float_t range = fGeomCam->GetMaxRadius()*1.05;
|
---|
1754 |
|
---|
1755 | if (!TestBit(kNoScale))
|
---|
1756 | {
|
---|
1757 | TArrow arr;
|
---|
1758 | arr.PaintArrow(-range*.9, -range*.9, -range*.6, -range*.9, 0.025);
|
---|
1759 | arr.PaintArrow(-range*.9, -range*.9, -range*.9, -range*.6, 0.025);
|
---|
1760 |
|
---|
1761 | TString text;
|
---|
1762 | text += (int)(range*.3);
|
---|
1763 | text += "mm";
|
---|
1764 |
|
---|
1765 | TText newtxt2;
|
---|
1766 | newtxt2.SetTextSize(0.04);
|
---|
1767 | newtxt2.PaintText(-range*.85, -range*.85, text);
|
---|
1768 |
|
---|
1769 | text = "";
|
---|
1770 | text += Form("%.2f", (float)((int)(range*.3*fGeomCam->GetConvMm2Deg()*10))/10);
|
---|
1771 | text += "\\circ";
|
---|
1772 | text = text.Strip(TString::kLeading);
|
---|
1773 |
|
---|
1774 | TLatex latex;
|
---|
1775 | latex.PaintLatex(-range*.85, -range*.75, 0, 0.04, text);
|
---|
1776 | }
|
---|
1777 |
|
---|
1778 | if (!TestBit(kNoLegend))
|
---|
1779 | {
|
---|
1780 | TPaveStats *stats = GetStatisticBox();
|
---|
1781 |
|
---|
1782 | const Float_t hndc = 0.92 - (stats ? stats->GetY1NDC() : 1);
|
---|
1783 | const Float_t H = (0.75-hndc)*range;
|
---|
1784 | const Float_t offset = hndc*range;
|
---|
1785 |
|
---|
1786 | const Float_t h = 2./kItemsLegend;
|
---|
1787 | const Float_t w = range/sqrt((float)(fNcells-2));
|
---|
1788 |
|
---|
1789 | TBox newbox;
|
---|
1790 | TText newtxt;
|
---|
1791 | newtxt.SetTextSize(0.03);
|
---|
1792 | newtxt.SetTextAlign(12);
|
---|
1793 | #if ROOT_VERSION_CODE > ROOT_VERSION(3,01,06)
|
---|
1794 | newtxt.SetBit(/*kNoContextMenu|*/kCannotPick);
|
---|
1795 | newbox.SetBit(/*kNoContextMenu|*/kCannotPick);
|
---|
1796 | #endif
|
---|
1797 |
|
---|
1798 | const Float_t step = (islog && min>0 ? log10(max/min) : max-min) / kItemsLegend;
|
---|
1799 | const Int_t firsts = step*3 < 1e-8 ? 8 : (Int_t)floor(log10(step*3));
|
---|
1800 | const TString opt = Form("%%.%if", firsts>0 ? 0 : TMath::Abs(firsts));
|
---|
1801 |
|
---|
1802 | for (Int_t i=0; i<kItemsLegend+1; i+=3)
|
---|
1803 | {
|
---|
1804 | Float_t val;
|
---|
1805 | if (islog && min>0)
|
---|
1806 | val = pow(10, step*i) * min;
|
---|
1807 | else
|
---|
1808 | val = min + step*i;
|
---|
1809 |
|
---|
1810 | //const bool dispexp = max-min>1.5 && fabs(val)>0.1 && fabs(val)<1e6;
|
---|
1811 | newtxt.PaintText(range+1.5*w, H*(i*h-1)-offset, Form(opt, val));
|
---|
1812 | }
|
---|
1813 |
|
---|
1814 | for (Int_t i=0; i<kItemsLegend; i++)
|
---|
1815 | {
|
---|
1816 | newbox.SetFillColor(gStyle->GetColorPalette(i));
|
---|
1817 | newbox.PaintBox(range, H*(i*h-1)-offset, range+w, H*((i+1)*h-1)-offset);
|
---|
1818 | }
|
---|
1819 | }
|
---|
1820 | }
|
---|
1821 |
|
---|
1822 | // ------------------------------------------------------------------------
|
---|
1823 | //
|
---|
1824 | // Save primitive as a C++ statement(s) on output stream out
|
---|
1825 | //
|
---|
1826 | void MHCamera::SavePrimitive(ostream &out, Option_t *opt)
|
---|
1827 | {
|
---|
1828 | gLog << err << "MHCamera::SavePrimitive: Must be rewritten!" << endl;
|
---|
1829 | /*
|
---|
1830 | if (!gROOT->ClassSaved(TCanvas::Class()))
|
---|
1831 | fDrawingPad->SavePrimitive(out, opt);
|
---|
1832 |
|
---|
1833 | out << " " << fDrawingPad->GetName() << "->SetWindowSize(";
|
---|
1834 | out << fDrawingPad->GetWw() << "," << fDrawingPad->GetWh() << ");" << endl;
|
---|
1835 | */
|
---|
1836 | }
|
---|
1837 |
|
---|
1838 | void MHCamera::SavePrimitive(ofstream &out, Option_t *)
|
---|
1839 | {
|
---|
1840 | MHCamera::SavePrimitive(static_cast<ostream&>(out), "");
|
---|
1841 | }
|
---|
1842 |
|
---|
1843 | // ------------------------------------------------------------------------
|
---|
1844 | //
|
---|
1845 | // compute the distance of a point (px,py) to the Camera
|
---|
1846 | // this functions needed for graphical primitives, that
|
---|
1847 | // means without this function you are not able to interact
|
---|
1848 | // with the graphical primitive with the mouse!!!
|
---|
1849 | //
|
---|
1850 | // All calcutations are done in pixel coordinates
|
---|
1851 | //
|
---|
1852 | Int_t MHCamera::DistancetoPrimitive(Int_t px, Int_t py)
|
---|
1853 | {
|
---|
1854 | if (fNcells<=1)
|
---|
1855 | return 999999;
|
---|
1856 |
|
---|
1857 | TPaveStats *box = (TPaveStats*)gPad->GetPrimitive("stats");
|
---|
1858 | if (box)
|
---|
1859 | {
|
---|
1860 | const Double_t w = box->GetY2NDC()-box->GetY1NDC();
|
---|
1861 | box->SetX1NDC(gStyle->GetStatX()-gStyle->GetStatW());
|
---|
1862 | box->SetY1NDC(gStyle->GetStatY()-w);
|
---|
1863 | box->SetX2NDC(gStyle->GetStatX());
|
---|
1864 | box->SetY2NDC(gStyle->GetStatY());
|
---|
1865 | }
|
---|
1866 |
|
---|
1867 | if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
|
---|
1868 | return TH1D::DistancetoPrimitive(px, py);
|
---|
1869 |
|
---|
1870 | const Bool_t issame = TString(GetDrawOption()).Contains("same", TString::kIgnoreCase);
|
---|
1871 |
|
---|
1872 | const Float_t maxr = (1-fGeomCam->GetConvMm2Deg())*fGeomCam->GetMaxRadius()/2;
|
---|
1873 | const Float_t conv = !issame ||
|
---|
1874 | gPad->GetX1()<-maxr || gPad->GetY1()<-maxr ||
|
---|
1875 | gPad->GetX2()> maxr || gPad->GetY2()>maxr ? 1 : fGeomCam->GetConvMm2Deg();
|
---|
1876 |
|
---|
1877 | if (GetPixelIndex(px, py, conv)>=0)
|
---|
1878 | return 0;
|
---|
1879 |
|
---|
1880 | if (!box)
|
---|
1881 | return 999999;
|
---|
1882 |
|
---|
1883 | const Int_t dist = box->DistancetoPrimitive(px, py);
|
---|
1884 | if (dist > TPad::GetMaxPickDistance())
|
---|
1885 | return 999999;
|
---|
1886 |
|
---|
1887 | gPad->SetSelected(box);
|
---|
1888 | return dist;
|
---|
1889 | }
|
---|
1890 |
|
---|
1891 | // ------------------------------------------------------------------------
|
---|
1892 | //
|
---|
1893 | //
|
---|
1894 | Int_t MHCamera::GetPixelIndex(Int_t px, Int_t py, Float_t conv) const
|
---|
1895 | {
|
---|
1896 | if (fNcells<=1)
|
---|
1897 | return -1;
|
---|
1898 |
|
---|
1899 | Int_t i;
|
---|
1900 | for (i=0; i<fNcells-2; i++)
|
---|
1901 | {
|
---|
1902 | MHexagon hex((*fGeomCam)[i]);
|
---|
1903 | if (hex.DistancetoPrimitive(px, py, conv)>0)
|
---|
1904 | continue;
|
---|
1905 |
|
---|
1906 | return i;
|
---|
1907 | }
|
---|
1908 | return -1;
|
---|
1909 | }
|
---|
1910 |
|
---|
1911 | // ------------------------------------------------------------------------
|
---|
1912 | //
|
---|
1913 | // Returns string containing info about the object at position (px,py).
|
---|
1914 | // Returned string will be re-used (lock in MT environment).
|
---|
1915 | //
|
---|
1916 | char *MHCamera::GetObjectInfo(Int_t px, Int_t py) const
|
---|
1917 | {
|
---|
1918 | if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
|
---|
1919 | return TH1D::GetObjectInfo(px, py);
|
---|
1920 |
|
---|
1921 | static char info[128];
|
---|
1922 |
|
---|
1923 | const Int_t idx=GetPixelIndex(px, py);
|
---|
1924 |
|
---|
1925 | if (idx<0)
|
---|
1926 | return TObject::GetObjectInfo(px, py);
|
---|
1927 |
|
---|
1928 | sprintf(info, "Software Pixel Idx: %d (Hardware Id=%d) c=%.1f <%s>",
|
---|
1929 | idx, idx+1, GetBinContent(idx+1), IsUsed(idx)?"on":"off");
|
---|
1930 | return info;
|
---|
1931 | }
|
---|
1932 |
|
---|
1933 | // ------------------------------------------------------------------------
|
---|
1934 | //
|
---|
1935 | // Add a MCamEvent which should be displayed when the user clicks on a
|
---|
1936 | // pixel.
|
---|
1937 | // Warning: The object MUST inherit from TObject AND MCamEvent
|
---|
1938 | //
|
---|
1939 | void MHCamera::AddNotify(TObject *obj)
|
---|
1940 | {
|
---|
1941 | // Make sure, that the object derives from MCamEvent!
|
---|
1942 | MCamEvent *evt = dynamic_cast<MCamEvent*>(obj);
|
---|
1943 | if (!evt)
|
---|
1944 | {
|
---|
1945 | gLog << err << "ERROR: MHCamera::AddNotify - TObject doesn't inherit from MCamEvent... ignored." << endl;
|
---|
1946 | return;
|
---|
1947 | }
|
---|
1948 |
|
---|
1949 | // Make sure, that it is deleted from the list too, if the obj is deleted
|
---|
1950 | obj->SetBit(kMustCleanup);
|
---|
1951 |
|
---|
1952 | // Add object to list
|
---|
1953 | fNotify->Add(obj);
|
---|
1954 | }
|
---|
1955 |
|
---|
1956 | // ------------------------------------------------------------------------
|
---|
1957 | //
|
---|
1958 | // Execute a mouse event on the camera
|
---|
1959 | //
|
---|
1960 | void MHCamera::ExecuteEvent(Int_t event, Int_t px, Int_t py)
|
---|
1961 | {
|
---|
1962 | if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
|
---|
1963 | {
|
---|
1964 | TH1D::ExecuteEvent(event, px, py);
|
---|
1965 | return;
|
---|
1966 | }
|
---|
1967 | //if (event==kMouseMotion && fStatusBar)
|
---|
1968 | // fStatusBar->SetText(GetObjectInfo(px, py), 0);
|
---|
1969 | if (event!=kButton1Down)
|
---|
1970 | return;
|
---|
1971 |
|
---|
1972 | const Int_t idx = GetPixelIndex(px, py);
|
---|
1973 | if (idx<0)
|
---|
1974 | return;
|
---|
1975 |
|
---|
1976 | gLog << all << GetTitle() << " <" << GetName() << ">" << dec << endl;
|
---|
1977 | gLog << "Software Pixel Idx: " << idx << endl;
|
---|
1978 | gLog << "Hardware Pixel Id: " << idx+1 << endl;
|
---|
1979 | gLog << "Contents: " << GetBinContent(idx+1);
|
---|
1980 | if (GetBinError(idx+1)>0)
|
---|
1981 | gLog << " +/- " << GetBinError(idx+1);
|
---|
1982 | gLog << " <" << (IsUsed(idx)?"on":"off") << "> n=" << fBinEntries[idx+1] << endl;
|
---|
1983 |
|
---|
1984 | if (fNotify && fNotify->GetSize()>0)
|
---|
1985 | {
|
---|
1986 | // FIXME: Is there a simpler and more convinient way?
|
---|
1987 |
|
---|
1988 | // The name which is created here depends on the instance of
|
---|
1989 | // MHCamera and on the pad on which it is drawn --> The name
|
---|
1990 | // is unique. For ExecuteEvent gPad is always correctly set.
|
---|
1991 | const TString name = Form("%p;%p;PixelContent", this, gPad);
|
---|
1992 |
|
---|
1993 | TCanvas *old = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(name);
|
---|
1994 | if (old)
|
---|
1995 | old->cd();
|
---|
1996 | else
|
---|
1997 | new TCanvas(name);
|
---|
1998 |
|
---|
1999 | /*
|
---|
2000 | TIter Next(gPad->GetListOfPrimitives());
|
---|
2001 | TObject *o;
|
---|
2002 | while (o=Next()) cout << o << ": " << o->GetName() << " " << o->IsA()->GetName() << endl;
|
---|
2003 | */
|
---|
2004 |
|
---|
2005 | // FIXME: Make sure, that the old histograms are really deleted.
|
---|
2006 | // Are they already deleted?
|
---|
2007 |
|
---|
2008 | // The dynamic_cast is necessary here: We cannot use ForEach
|
---|
2009 | TIter Next(fNotify);
|
---|
2010 | MCamEvent *evt;
|
---|
2011 | while ((evt=dynamic_cast<MCamEvent*>(Next())))
|
---|
2012 | evt->DrawPixelContent(idx);
|
---|
2013 |
|
---|
2014 | gPad->Modified();
|
---|
2015 | gPad->Update();
|
---|
2016 | }
|
---|
2017 | }
|
---|
2018 |
|
---|
2019 | UInt_t MHCamera::GetNumPixels() const
|
---|
2020 | {
|
---|
2021 | return fGeomCam ? fGeomCam->GetNumPixels() : 0;
|
---|
2022 | }
|
---|
2023 |
|
---|
2024 | TH1 *MHCamera::DrawCopy() const
|
---|
2025 | {
|
---|
2026 | gPad=NULL;
|
---|
2027 | return TH1D::DrawCopy(fName+";cpy");
|
---|
2028 | }
|
---|
2029 |
|
---|
2030 | // --------------------------------------------------------------------------
|
---|
2031 | //
|
---|
2032 | // Draw a projection of MHCamera onto the y-axis values. Depending on the
|
---|
2033 | // variable fit, the following fits are performed:
|
---|
2034 | //
|
---|
2035 | // 0: No fit, simply draw the projection
|
---|
2036 | // 1: Single Gauss (for distributions flat-fielded over the whole camera)
|
---|
2037 | // 2: Double Gauss (for distributions different for inner and outer pixels)
|
---|
2038 | // 3: Triple Gauss (for distributions with inner, outer pixels and outliers)
|
---|
2039 | // 4: flat (for the probability distributions)
|
---|
2040 | // (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are
|
---|
2041 | // drawn separately, for inner and outer pixels.
|
---|
2042 | // 5: Fit Inner and Outer pixels separately by a single Gaussian
|
---|
2043 | // (only for MAGIC cameras)
|
---|
2044 | // 6: Fit Inner and Outer pixels separately by a single Gaussian and display
|
---|
2045 | // additionally the two camera halfs separately (for MAGIC camera)
|
---|
2046 | // 7: Single Gauss with TLegend to show the meaning of the colours
|
---|
2047 | //
|
---|
2048 | void MHCamera::DrawProjection(Int_t fit) const
|
---|
2049 | {
|
---|
2050 | TArrayI inner(1);
|
---|
2051 | inner[0] = 0;
|
---|
2052 |
|
---|
2053 | TArrayI outer(1);
|
---|
2054 | outer[0] = 1;
|
---|
2055 |
|
---|
2056 | if (fit==5 || fit==6)
|
---|
2057 | {
|
---|
2058 | if (GetGeomCam().InheritsFrom("MGeomCamMagic"))
|
---|
2059 | {
|
---|
2060 | TArrayI s0(6);
|
---|
2061 | s0[0] = 6;
|
---|
2062 | s0[1] = 1;
|
---|
2063 | s0[2] = 2;
|
---|
2064 | s0[3] = 3;
|
---|
2065 | s0[4] = 4;
|
---|
2066 | s0[5] = 5;
|
---|
2067 |
|
---|
2068 | TArrayI s1(3);
|
---|
2069 | s1[0] = 6;
|
---|
2070 | s1[1] = 1;
|
---|
2071 | s1[2] = 2;
|
---|
2072 |
|
---|
2073 | TArrayI s2(3);
|
---|
2074 | s2[0] = 3;
|
---|
2075 | s2[1] = 4;
|
---|
2076 | s2[2] = 5;
|
---|
2077 |
|
---|
2078 | gPad->Clear();
|
---|
2079 | TVirtualPad *pad = gPad;
|
---|
2080 | pad->Divide(2,1);
|
---|
2081 |
|
---|
2082 | TH1D *inout[2];
|
---|
2083 | inout[0] = ProjectionS(s0, inner, "Inner");
|
---|
2084 | inout[1] = ProjectionS(s0, outer, "Outer");
|
---|
2085 |
|
---|
2086 | inout[0]->SetDirectory(NULL);
|
---|
2087 | inout[1]->SetDirectory(NULL);
|
---|
2088 |
|
---|
2089 | for (int i=0; i<2; i++)
|
---|
2090 | {
|
---|
2091 | pad->cd(i+1);
|
---|
2092 | gPad->SetBorderMode(0);
|
---|
2093 |
|
---|
2094 | inout[i]->SetLineColor(kRed+i);
|
---|
2095 | inout[i]->SetBit(kCanDelete);
|
---|
2096 | inout[i]->Draw();
|
---|
2097 | inout[i]->Fit("gaus","Q");
|
---|
2098 |
|
---|
2099 | if (fit == 6)
|
---|
2100 | {
|
---|
2101 | TH1D *half[2];
|
---|
2102 | half[0] = ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
|
---|
2103 | half[1] = ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
|
---|
2104 |
|
---|
2105 | for (int j=0; j<2; j++)
|
---|
2106 | {
|
---|
2107 | half[j]->SetLineColor(kRed+i+2*j+1);
|
---|
2108 | half[j]->SetDirectory(NULL);
|
---|
2109 | half[j]->SetBit(kCanDelete);
|
---|
2110 | half[j]->Draw("same");
|
---|
2111 | }
|
---|
2112 | }
|
---|
2113 |
|
---|
2114 | }
|
---|
2115 | }
|
---|
2116 | return;
|
---|
2117 | }
|
---|
2118 |
|
---|
2119 | TH1D *obj2 = (TH1D*)Projection(GetName());
|
---|
2120 | obj2->SetDirectory(0);
|
---|
2121 | obj2->Draw();
|
---|
2122 | obj2->SetBit(kCanDelete);
|
---|
2123 |
|
---|
2124 | if (fit == 0)
|
---|
2125 | return;
|
---|
2126 |
|
---|
2127 | if (GetGeomCam().InheritsFrom("MGeomCamMagic"))
|
---|
2128 | {
|
---|
2129 | TArrayI s0(3);
|
---|
2130 | s0[0] = 6;
|
---|
2131 | s0[1] = 1;
|
---|
2132 | s0[2] = 2;
|
---|
2133 |
|
---|
2134 | TArrayI s1(3);
|
---|
2135 | s1[0] = 3;
|
---|
2136 | s1[1] = 4;
|
---|
2137 | s1[2] = 5;
|
---|
2138 |
|
---|
2139 | TH1D *halfInOut[4];
|
---|
2140 |
|
---|
2141 | // Just to get the right (maximum) binning
|
---|
2142 | halfInOut[0] = ProjectionS(s0, inner, "Sector 6-1-2 Inner");
|
---|
2143 | halfInOut[1] = ProjectionS(s1, inner, "Sector 3-4-5 Inner");
|
---|
2144 | halfInOut[2] = ProjectionS(s0, outer, "Sector 6-1-2 Outer");
|
---|
2145 | halfInOut[3] = ProjectionS(s1, outer, "Sector 3-4-5 Outer");
|
---|
2146 |
|
---|
2147 | TLegend *leg = new TLegend(0.05,0.65,0.35,0.9);
|
---|
2148 |
|
---|
2149 | for (int i=0; i<4; i++)
|
---|
2150 | {
|
---|
2151 | halfInOut[i]->SetLineColor(kRed+i);
|
---|
2152 | halfInOut[i]->SetDirectory(0);
|
---|
2153 | halfInOut[i]->SetBit(kCanDelete);
|
---|
2154 | halfInOut[i]->Draw("same");
|
---|
2155 | leg->AddEntry(halfInOut[i],halfInOut[i]->GetTitle(),"l");
|
---|
2156 | }
|
---|
2157 |
|
---|
2158 | if (fit==7)
|
---|
2159 | leg->Draw();
|
---|
2160 |
|
---|
2161 | gPad->Modified();
|
---|
2162 | gPad->Update();
|
---|
2163 | }
|
---|
2164 |
|
---|
2165 | const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
|
---|
2166 | const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
|
---|
2167 | const Double_t integ = obj2->Integral("width")/2.5;
|
---|
2168 | const Double_t mean = obj2->GetMean();
|
---|
2169 | const Double_t rms = obj2->GetRMS();
|
---|
2170 | const Double_t width = max-min;
|
---|
2171 |
|
---|
2172 | const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
|
---|
2173 | "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
|
---|
2174 |
|
---|
2175 | const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
|
---|
2176 | "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
|
---|
2177 | "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
|
---|
2178 |
|
---|
2179 | TF1 *f=0;
|
---|
2180 | switch (fit)
|
---|
2181 | {
|
---|
2182 | case 1:
|
---|
2183 | f = new TF1("sgaus", "gaus(0)", min, max);
|
---|
2184 | f->SetLineColor(kYellow);
|
---|
2185 | f->SetBit(kCanDelete);
|
---|
2186 | f->SetParNames("Area", "#mu", "#sigma");
|
---|
2187 | f->SetParameters(integ/rms, mean, rms);
|
---|
2188 | f->SetParLimits(0, 0, integ);
|
---|
2189 | f->SetParLimits(1, min, max);
|
---|
2190 | f->SetParLimits(2, 0, width/1.5);
|
---|
2191 |
|
---|
2192 | obj2->Fit(f, "QLR");
|
---|
2193 | break;
|
---|
2194 |
|
---|
2195 | case 2:
|
---|
2196 | f = new TF1("dgaus",dgausformula.Data(),min,max);
|
---|
2197 | f->SetLineColor(kYellow);
|
---|
2198 | f->SetBit(kCanDelete);
|
---|
2199 | f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
|
---|
2200 | f->SetParameters(integ,(min+mean)/2.,width/4.,
|
---|
2201 | integ/width/2.,(max+mean)/2.,width/4.);
|
---|
2202 | // The left-sided Gauss
|
---|
2203 | f->SetParLimits(0,integ-1.5 , integ+1.5);
|
---|
2204 | f->SetParLimits(1,min+(width/10.), mean);
|
---|
2205 | f->SetParLimits(2,0 , width/2.);
|
---|
2206 | // The right-sided Gauss
|
---|
2207 | f->SetParLimits(3,0 , integ);
|
---|
2208 | f->SetParLimits(4,mean, max-(width/10.));
|
---|
2209 | f->SetParLimits(5,0 , width/2.);
|
---|
2210 | obj2->Fit(f,"QLRM");
|
---|
2211 | break;
|
---|
2212 |
|
---|
2213 | case 3:
|
---|
2214 | f = new TF1("tgaus",tgausformula.Data(),min,max);
|
---|
2215 | f->SetLineColor(kYellow);
|
---|
2216 | f->SetBit(kCanDelete);
|
---|
2217 | f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
|
---|
2218 | "A_{2}","#mu_{2}","#sigma_{2}",
|
---|
2219 | "A_{3}","#mu_{3}","#sigma_{3}");
|
---|
2220 | f->SetParameters(integ,(min+mean)/2,width/4.,
|
---|
2221 | integ/width/3.,(max+mean)/2.,width/4.,
|
---|
2222 | integ/width/3.,mean,width/2.);
|
---|
2223 | // The left-sided Gauss
|
---|
2224 | f->SetParLimits(0,integ-1.5,integ+1.5);
|
---|
2225 | f->SetParLimits(1,min+(width/10.),mean);
|
---|
2226 | f->SetParLimits(2,width/15.,width/2.);
|
---|
2227 | // The right-sided Gauss
|
---|
2228 | f->SetParLimits(3,0.,integ);
|
---|
2229 | f->SetParLimits(4,mean,max-(width/10.));
|
---|
2230 | f->SetParLimits(5,width/15.,width/2.);
|
---|
2231 | // The Gauss describing the outliers
|
---|
2232 | f->SetParLimits(6,0.,integ);
|
---|
2233 | f->SetParLimits(7,min,max);
|
---|
2234 | f->SetParLimits(8,width/4.,width/1.5);
|
---|
2235 | obj2->Fit(f,"QLRM");
|
---|
2236 | break;
|
---|
2237 |
|
---|
2238 | case 4:
|
---|
2239 | obj2->Fit("pol0", "Q");
|
---|
2240 | obj2->GetFunction("pol0")->SetLineColor(kYellow);
|
---|
2241 | break;
|
---|
2242 |
|
---|
2243 | case 9:
|
---|
2244 | break;
|
---|
2245 |
|
---|
2246 | default:
|
---|
2247 | obj2->Fit("gaus", "Q");
|
---|
2248 | obj2->GetFunction("gaus")->SetLineColor(kYellow);
|
---|
2249 | break;
|
---|
2250 | }
|
---|
2251 | }
|
---|
2252 |
|
---|
2253 | // --------------------------------------------------------------------------
|
---|
2254 | //
|
---|
2255 | // Draw a projection of MHCamera vs. the radius from the central pixel.
|
---|
2256 | //
|
---|
2257 | // The inner and outer pixels are drawn separately, both fitted by a polynomial
|
---|
2258 | // of grade 1.
|
---|
2259 | //
|
---|
2260 | void MHCamera::DrawRadialProfile() const
|
---|
2261 | {
|
---|
2262 | TProfile *obj2 = (TProfile*)RadialProfile(GetName());
|
---|
2263 | obj2->SetDirectory(0);
|
---|
2264 | obj2->Draw();
|
---|
2265 | obj2->SetBit(kCanDelete);
|
---|
2266 |
|
---|
2267 | if (GetGeomCam().InheritsFrom("MGeomCamMagic"))
|
---|
2268 | {
|
---|
2269 | TArrayI s0(6);
|
---|
2270 | s0[0] = 1;
|
---|
2271 | s0[1] = 2;
|
---|
2272 | s0[2] = 3;
|
---|
2273 | s0[3] = 4;
|
---|
2274 | s0[4] = 5;
|
---|
2275 | s0[5] = 6;
|
---|
2276 |
|
---|
2277 | TArrayI inner(1);
|
---|
2278 | inner[0] = 0;
|
---|
2279 |
|
---|
2280 | TArrayI outer(1);
|
---|
2281 | outer[0] = 1;
|
---|
2282 |
|
---|
2283 | // Just to get the right (maximum) binning
|
---|
2284 | TProfile *half[2];
|
---|
2285 | half[0] = RadialProfileS(s0, inner,Form("%sInner",GetName()));
|
---|
2286 | half[1] = RadialProfileS(s0, outer,Form("%sOuter",GetName()));
|
---|
2287 |
|
---|
2288 | for (Int_t i=0; i<2; i++)
|
---|
2289 | {
|
---|
2290 | Double_t min = GetGeomCam().GetMinRadius(i);
|
---|
2291 | Double_t max = GetGeomCam().GetMaxRadius(i);
|
---|
2292 |
|
---|
2293 | half[i]->SetLineColor(kRed+i);
|
---|
2294 | half[i]->SetDirectory(0);
|
---|
2295 | half[i]->SetBit(kCanDelete);
|
---|
2296 | half[i]->Draw("same");
|
---|
2297 | half[i]->Fit("pol1","Q","",min,max);
|
---|
2298 | half[i]->GetFunction("pol1")->SetLineColor(kRed+i);
|
---|
2299 | half[i]->GetFunction("pol1")->SetLineWidth(1);
|
---|
2300 | }
|
---|
2301 | }
|
---|
2302 | }
|
---|
2303 |
|
---|
2304 | // --------------------------------------------------------------------------
|
---|
2305 | //
|
---|
2306 | // Draw a projection of MHCamera vs. the azimuth angle inside the camera.
|
---|
2307 | //
|
---|
2308 | // The inner and outer pixels are drawn separately.
|
---|
2309 | // The general azimuth profile is fitted by a straight line
|
---|
2310 | //
|
---|
2311 | void MHCamera::DrawAzimuthProfile() const
|
---|
2312 | {
|
---|
2313 | TProfile *obj2 = (TProfile*)AzimuthProfile(GetName());
|
---|
2314 | obj2->SetDirectory(0);
|
---|
2315 | obj2->Draw();
|
---|
2316 | obj2->SetBit(kCanDelete);
|
---|
2317 | obj2->Fit("pol0","Q","");
|
---|
2318 | obj2->GetFunction("pol0")->SetLineWidth(1);
|
---|
2319 |
|
---|
2320 | if (GetGeomCam().InheritsFrom("MGeomCamMagic"))
|
---|
2321 | {
|
---|
2322 | TArrayI inner(1);
|
---|
2323 | inner[0] = 0;
|
---|
2324 |
|
---|
2325 | TArrayI outer(1);
|
---|
2326 | outer[0] = 1;
|
---|
2327 |
|
---|
2328 | // Just to get the right (maximum) binning
|
---|
2329 | TProfile *half[2];
|
---|
2330 | half[0] = AzimuthProfileA(inner,Form("%sInner",GetName()));
|
---|
2331 | half[1] = AzimuthProfileA(outer,Form("%sOuter",GetName()));
|
---|
2332 |
|
---|
2333 | for (Int_t i=0; i<2; i++)
|
---|
2334 | {
|
---|
2335 | half[i]->SetLineColor(kRed+i);
|
---|
2336 | half[i]->SetDirectory(0);
|
---|
2337 | half[i]->SetBit(kCanDelete);
|
---|
2338 | half[i]->SetMarkerSize(0.5);
|
---|
2339 | half[i]->Draw("same");
|
---|
2340 | }
|
---|
2341 | }
|
---|
2342 | }
|
---|
2343 |
|
---|
2344 | // --------------------------------------------------------------------------
|
---|
2345 | //
|
---|
2346 | // Draw the MHCamera into the MStatusDisplay:
|
---|
2347 | //
|
---|
2348 | // 1) Draw it as histogram (MHCamera::DrawCopy("hist")
|
---|
2349 | // 2) Draw it as a camera, with MHCamera::SetPrettyPalette() set.
|
---|
2350 | // 3) If "rad" is not zero, draw its values vs. the radius from the camera center.
|
---|
2351 | // (DrawRadialProfile())
|
---|
2352 | // 4) Depending on the variable "fit", draw the values projection on the y-axis
|
---|
2353 | // (DrawProjection()):
|
---|
2354 | // 0: don't draw
|
---|
2355 | // 1: Draw fit to Single Gauss (for distributions flat-fielded over the whole camera)
|
---|
2356 | // 2: Draw and fit to Double Gauss (for distributions different for inner and outer pixels)
|
---|
2357 | // 3: Draw and fit to Triple Gauss (for distributions with inner, outer pixels and outliers)
|
---|
2358 | // 4: Draw and fit to Polynomial grade 0: (for the probability distributions)
|
---|
2359 | // >4: Draw and don;t fit.
|
---|
2360 | //
|
---|
2361 | void MHCamera::CamDraw(TCanvas &c, const Int_t x, const Int_t y,
|
---|
2362 | const Int_t fit, const Int_t rad, const Int_t azi,
|
---|
2363 | TObject *notify)
|
---|
2364 | {
|
---|
2365 | c.cd(x);
|
---|
2366 | gPad->SetBorderMode(0);
|
---|
2367 | gPad->SetTicks();
|
---|
2368 | MHCamera *obj1=(MHCamera*)DrawCopy("hist");
|
---|
2369 | obj1->SetDirectory(NULL);
|
---|
2370 |
|
---|
2371 | if (notify)
|
---|
2372 | obj1->AddNotify(notify);
|
---|
2373 |
|
---|
2374 | c.cd(x+y);
|
---|
2375 | gPad->SetBorderMode(0);
|
---|
2376 | obj1->SetPrettyPalette();
|
---|
2377 | obj1->Draw();
|
---|
2378 |
|
---|
2379 | Int_t cnt = 2;
|
---|
2380 |
|
---|
2381 | if (rad)
|
---|
2382 | {
|
---|
2383 | c.cd(x+2*y);
|
---|
2384 | gPad->SetBorderMode(0);
|
---|
2385 | gPad->SetTicks();
|
---|
2386 | DrawRadialProfile();
|
---|
2387 | cnt++;
|
---|
2388 | }
|
---|
2389 |
|
---|
2390 | if (azi)
|
---|
2391 | {
|
---|
2392 | c.cd(x+cnt*y);
|
---|
2393 | gPad->SetBorderMode(0);
|
---|
2394 | gPad->SetTicks();
|
---|
2395 | DrawAzimuthProfile();
|
---|
2396 | cnt++;
|
---|
2397 | }
|
---|
2398 |
|
---|
2399 | if (!fit)
|
---|
2400 | return;
|
---|
2401 |
|
---|
2402 | c.cd(x + cnt*y);
|
---|
2403 | gPad->SetBorderMode(0);
|
---|
2404 | gPad->SetTicks();
|
---|
2405 | DrawProjection(fit);
|
---|
2406 | }
|
---|