source: trunk/MagicSoft/Mars/mhist/MHCamera.cc@ 4883

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