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

Last change on this file since 8281 was 8281, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 68.9 KB
Line 
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
92ClassImp(MHCamera);
93
94using namespace std;
95
96void 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//
134MHCamera::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//
145MHCamera::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
159void 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//
187MHCamera::~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//
200Bool_t MHCamera::MatchSector(Int_t idx, const TArrayI &sector, 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// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
222Int_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// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
259Int_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//
285Int_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//
307void 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//
325Stat_t MHCamera::GetMeanSectors(const TArrayI &sector, 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//
354Stat_t MHCamera::GetMedianSectors(const TArrayI &sector, 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//
383Stat_t MHCamera::GetRmsSectors(const TArrayI &sector, 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//
424Stat_t MHCamera::GetDevSectors(const TArrayI &sector, 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//
453Double_t MHCamera::GetMinimumSectors(const TArrayI &sector, 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//
482Double_t MHCamera::GetMaximumSectors(const TArrayI &sector, 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//
541void 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//
582TObject *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//
646TH1D *MHCamera::ProjectionS(const TArrayI &sector, 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//
741TProfile *MHCamera::RadialProfileS(const TArrayI &sector, 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//
827TProfile *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//
890void 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//
932void 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//
998void 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//
1014void 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//
1036void 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
1140void 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//
1167void 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//
1200void 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//
1219void 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//
1239void 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//
1261void 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//
1302void 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//
1324void 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
1341Stat_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//
1361Stat_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//
1400void 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//
1470void 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//
1492void 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//
1518void 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//
1552void 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//
1585void 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//
1612void 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//
1632void 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//
1657void 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//
1681void 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//
1706Int_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
1735TPaveStats *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//
1751void 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//
1826void 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
1838void 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//
1852Int_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//
1894Int_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//
1916char *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//
1939void 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//
1960void 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
2019UInt_t MHCamera::GetNumPixels() const
2020{
2021 return fGeomCam ? fGeomCam->GetNumPixels() : 0;
2022}
2023
2024TH1 *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//
2048void 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//
2260void 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//
2311void 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//
2361void 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}
Note: See TracBrowser for help on using the repository browser.