source: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc@ 9029

Last change on this file since 9029 was 8996, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 31.6 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, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHAlpha
28//
29// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
30// check the result when it is filles in the MStatusDisplay
31// For more information see MHFalseSource::FitSignificance
32//
33// For convinience (fit) the output significance is stored in a
34// container in the parlisrt
35//
36// PRELIMINARY!
37//
38//
39// ToDo:
40// =====
41//
42// - Replace time-binning by histogram (which is enhanced bin-wise)
43//
44//
45//////////////////////////////////////////////////////////////////////////////
46#include "MHAlpha.h"
47
48#include <TF1.h>
49#include <TGraph.h>
50#include <TStyle.h>
51#include <TLatex.h>
52#include <TCanvas.h>
53#include <TPaveStats.h>
54#include <TStopwatch.h>
55
56#include "MSrcPosCam.h"
57#include "MHillas.h"
58#include "MHillasSrc.h"
59#include "MTime.h"
60#include "MObservatory.h"
61#include "MPointingPos.h"
62#include "MAstroSky2Local.h"
63#include "MStatusDisplay.h"
64#include "MParameters.h"
65#include "MHMatrix.h"
66
67#include "MBinning.h"
68#include "MParList.h"
69
70#include "MLog.h"
71#include "MLogManip.h"
72
73ClassImp(MHAlpha);
74
75using namespace std;
76
77// --------------------------------------------------------------------------
78//
79// Default Constructor
80//
81MHAlpha::MHAlpha(const char *name, const char *title)
82 : fNameParameter("MHillasSrc"), fParameter(0),
83 fOffData(0), fResult(0), fEnergy(0), fBin(0),
84 fPointPos(0), fTimeEffOn(0), fTime(0), fNumTimeBins(10),
85 fHillas(0), fMatrix(0), fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE),
86 fSkipHistEnergy(kFALSE), fForceUsingSize(kFALSE)
87{
88 //
89 // set the name and title of this object
90 //
91 fName = name ? name : "MHAlpha";
92 fTitle = title ? title : "Alpha plot";
93
94 fHist.SetName("Alpha");
95 fHist.SetTitle("Alpha");
96 fHist.SetXTitle("\\Theta [deg]");
97 //fHist.SetYTitle("E_{est} [GeV]");
98 fHist.SetZTitle("|\\alpha| [\\circ]");
99 fHist.SetDirectory(NULL);
100 fHist.UseCurrentStyle();
101
102 // Main histogram
103 fHistTime.SetName("Alpha");
104 fHistTime.SetXTitle("|\\alpha| [\\circ]");
105 fHistTime.SetYTitle("Counts");
106 fHistTime.UseCurrentStyle();
107 fHistTime.SetDirectory(NULL);
108
109
110 fHEnergy.SetName("Excess");
111 //fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
112 //fHEnergy.SetXTitle("E_{est} [GeV]");
113 fHEnergy.SetYTitle("N_{exc}");
114 fHEnergy.SetDirectory(NULL);
115 fHEnergy.UseCurrentStyle();
116
117 fHTheta.SetName("ExcessTheta");
118 fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
119 fHTheta.SetXTitle("\\Theta [\\circ]");
120 fHTheta.SetYTitle("N_{exc}");
121 fHTheta.SetDirectory(NULL);
122 fHTheta.UseCurrentStyle();
123 fHTheta.SetMinimum(0);
124
125 // effective on time versus time
126 fHTime.SetName("ExcessTime");
127 fHTime.SetTitle(" N_{exc} vs. Time ");
128 fHTime.SetXTitle("Time");
129 fHTime.SetYTitle("N_{exc} [s]");
130 fHTime.UseCurrentStyle();
131 fHTime.SetDirectory(NULL);
132 fHTime.GetYaxis()->SetTitleOffset(1.2);
133 fHTime.GetXaxis()->SetLabelSize(0.033);
134 fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
135 fHTime.GetXaxis()->SetTimeDisplay(1);
136 fHTime.SetMinimum(0);
137 fHTime.Sumw2();
138
139 MBinning binsa, binse, binst;
140 binsa.SetEdges(18, 0, 90);
141 binse.SetEdgesLog(15, 10, 100000);
142 binst.SetEdgesASin(67, -0.005, 0.665);
143 binse.Apply(fHEnergy);
144 binst.Apply(fHTheta);
145 binsa.Apply(fHistTime);
146
147 MH::SetBinning(&fHist, &binst, &binse, &binsa);
148}
149
150Float_t MHAlpha::FitEnergyBins(Bool_t paint)
151{
152 // Do not store the 'final' result in fFit
153 MAlphaFitter fit(fFit);
154
155 const Int_t n = fHist.GetNbinsY();
156
157 fHEnergy.SetEntries(0);
158
159 Float_t mean = 0;
160 Int_t num = 0;
161 for (int i=1; i<=n; i++)
162 {
163 if (!fit.FitEnergy(fHist, fOffData, i))
164 continue;
165
166 // FIXME: Calculate UL!
167 if (fit.GetEventsExcess()<=0)
168 continue;
169
170 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
171 fHEnergy.SetBinError(i, fit.GetErrorExcess());
172
173 mean += fit.GetEventsExcess()*fit.GetEventsExcess()/fit.GetErrorExcess()/fit.GetErrorExcess();
174 num++;
175 }
176 return TMath::Sqrt(mean)/num;
177}
178
179void MHAlpha::FitThetaBins(Bool_t paint)
180{
181 // Do not store the 'final' result in fFit
182 MAlphaFitter fit(fFit);
183
184 const Int_t n = fHist.GetNbinsX();
185
186 fHTheta.SetEntries(0);
187
188 for (int i=1; i<=n; i++)
189 {
190 if (!fit.FitTheta(fHist, fOffData, i))
191 continue;
192
193 // FIXME: Calculate UL!
194 if (fit.GetEventsExcess()<=0)
195 continue;
196
197 fHTheta.SetBinContent(i, fit.GetEventsExcess());
198 fHTheta.SetBinError(i, fit.GetErrorExcess());
199 }
200}
201
202// --------------------------------------------------------------------------
203//
204// Return the value from fParemeter which should be filled into the plots
205//
206Double_t MHAlpha::GetVal() const
207{
208 return static_cast<const MHillasSrc*>(fParameter)->GetAlpha();
209}
210
211
212// --------------------------------------------------------------------------
213//
214// Store the pointer to the parameter container storing the plotted value
215// (here MHillasSrc) in fParameter.
216//
217// return whether it was found or not.
218//
219Bool_t MHAlpha::GetParameter(const MParList &pl)
220{
221 fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MHillasSrc");
222 if (fParameter)
223 return kTRUE;
224
225 *fLog << err << fNameParameter << " [MHillasSrc] not found... abort." << endl;
226 return kFALSE;
227}
228
229Bool_t MHAlpha::SetupFill(const MParList *pl)
230{
231 fHist.Reset();
232 fHEnergy.Reset();
233 fHTheta.Reset();
234 fHTime.Reset();
235
236 const TString off(Form("%sOff", fName.Data()));
237 if (fName!=off && fOffData==NULL)
238 {
239 const TString desc(Form("%s [%s] found... using ", off.Data(), ClassName()));
240 MHAlpha *hoff = (MHAlpha*)pl->FindObject(off, ClassName());
241 if (!hoff)
242 *fLog << inf << "No " << desc << "current data only!" << endl;
243 else
244 {
245 *fLog << inf << desc << "on-off mode!" << endl;
246 SetOffData(*hoff);
247 }
248 }
249
250 if (!GetParameter(*pl))
251 return kFALSE;
252
253 fHillas = 0;
254 fEnergy = fForceUsingSize ? 0 : (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD");
255 if (!fEnergy)
256 {
257 *fLog << warn << "MEnergyEst [MParameterD] not found... " << flush;
258
259 if (!fHillas)
260 fHillas = (MHillas*)pl->FindObject("MHillas");
261 if (fHillas)
262 *fLog << "using SIZE instead!" << endl;
263 else
264 *fLog << "ignored." << endl;
265
266 fHEnergy.SetTitle(" N_{exc} vs. Size ");
267 fHEnergy.SetXTitle("Size [phe]");
268 fHist.SetYTitle("Size [phe]");
269 }
270 else
271 {
272 fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
273 fHEnergy.SetXTitle("E_{est} [GeV]");
274 fHist.SetYTitle("E_{est} [GeV]");
275 }
276
277 fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
278 if (!fPointPos)
279 *fLog << warn << "MPointingPos not found... ignored." << endl;
280
281 fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
282 if (!fTimeEffOn)
283 *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
284 else
285 *fTimeEffOn = MTime(); // FIXME: How to do it different?
286
287 fTime = (MTime*)pl->FindObject("MTime");
288 if (!fTime)
289 *fLog << warn << "MTime not found... ignored." << endl;
290
291 fResult = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MinimizationValue");
292 if (!fResult)
293 return kFALSE;
294 fSigma = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "GaussSigma");
295 if (!fSigma)
296 return kFALSE;
297 fBin = (MParameterI*)const_cast<MParList*>(pl)->FindCreateObj("MParameterI", "Bin");
298 if (!fBin)
299 return kFALSE;
300
301 //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
302 //if (!fExcess)
303 // return kFALSE;
304
305 fLastTime = MTime();
306 fNumRebin = fNumTimeBins;
307
308 MBinning binst, binse, binsa;
309 binst.SetEdges(fOffData ? *fOffData : fHist, 'x');
310 binse.SetEdges(fOffData ? *fOffData : fHist, 'y');
311 binsa.SetEdges(fOffData ? *fOffData : fHist, 'z');
312 if (!fOffData)
313 {
314 if (!fPointPos)
315 binst.SetEdges(1, 0, 60);
316 else
317 binst.SetEdges(*pl, "BinningTheta");
318
319 if (!fEnergy && !fHillas)
320 binse.SetEdges(1, 10, 100000);
321 else
322 if (fEnergy)
323 binse.SetEdges(*pl, "BinningEnergyEst");
324 else
325 binse.SetEdges(*pl, "BinningSize");
326
327 binsa.SetEdges(*pl, Form("Binning%s", ClassName()+2));
328 }
329 else
330 {
331 fHEnergy.SetTitle(fOffData->GetTitle());
332 fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle());
333 fHist.SetYTitle(fOffData->GetYaxis()->GetTitle());
334 }
335
336 binse.Apply(fHEnergy);
337 binst.Apply(fHTheta);
338 binsa.Apply(fHistTime);
339 MH::SetBinning(&fHist, &binst, &binse, &binsa);
340
341 MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
342 if (!fit)
343 *fLog << warn << "MAlphaFitter not found... ignored." << endl;
344 else
345 fit->Copy(fFit);
346
347 *fLog << inf;
348 fFit.Print();
349
350 return kTRUE;
351}
352
353void MHAlpha::InitAlphaTime(const MTime &t)
354{
355 //
356 // If this is the first call we have to initialize the time-histogram
357 //
358 MBinning bins;
359 bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
360 bins.Apply(fHTime);
361
362 fLastTime=t;
363}
364
365void MHAlpha::UpdateAlphaTime(Bool_t final)
366{
367 if (!fTimeEffOn)
368 return;
369
370 if (!final)
371 {
372 if (fLastTime==MTime() && *fTimeEffOn!=MTime())
373 {
374 InitAlphaTime(*fTimeEffOn);
375 return;
376 }
377
378 if (fLastTime!=*fTimeEffOn)
379 {
380 fLastTime=*fTimeEffOn;
381 fNumRebin--;
382 }
383
384 if (fNumRebin>0)
385 return;
386 }
387
388 // Check new 'last time'
389 MTime *time = final ? fTime : fTimeEffOn;
390
391 if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
392 {
393 *fLog << warn << "WARNING - New time-stamp " << *time << " lower" << endl;
394 *fLog << "than upper edge of histogram... skipped." << endl;
395 *fLog << "This should not happen. Maybe you started you eventloop with" << endl;
396 *fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
397 fNumRebin++;
398 return;
399 }
400
401 // Fit time histogram
402 MAlphaFitter fit(fFit);
403
404 TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E") : 0;
405 const Bool_t rc = fit.ScaleAndFit(fHistTime, h);
406
407 if (h)
408 delete h;
409
410 if (!rc)
411 return;
412
413 // Reset Histogram
414 fHistTime.Reset();
415 fHistTime.SetEntries(0);
416
417 //
418 // Prepare Histogram
419 //
420 if (final)
421 time->Plus1ns();
422
423 // Enhance binning
424 MBinning bins;
425 bins.SetEdges(fHTime, 'x');
426 bins.AddEdge(time->GetAxisTime());
427 bins.Apply(fHTime);
428
429 //
430 // Fill histogram
431 //
432 // Get number of bins
433 const Int_t n = fHTime.GetNbinsX();
434
435 fHTime.SetBinContent(n, fit.GetEventsExcess());
436 fHTime.SetBinError(n, fit.GetErrorExcess());
437
438 *fLog << all << *fTimeEffOn << " (" << n << "): " << fit.GetEventsExcess() << endl;
439
440 fNumRebin = fNumTimeBins;
441}
442
443void MHAlpha::SetBin(Int_t ibin)
444{
445 // Is this necessary?
446 // Could be speed up up searching for it only once.
447 const Float_t max = fFit.GetSignalIntegralMax();
448 const Int_t bin0 = fHist.GetZaxis()->FindFixBin(max);
449
450 const Int_t nbinsx = fHist.GetNbinsX();
451 const Int_t nbinsy = fHist.GetNbinsY();
452 const Int_t nxy = (nbinsx+2)*(nbinsy+2);
453
454 const Int_t binz = ibin/nxy;
455
456 const Bool_t issignal = binz>0 && binz<bin0;
457
458 fBin->SetVal(issignal ? binz : -binz);
459}
460
461// --------------------------------------------------------------------------
462//
463// Fill the histogram. For details see the code or the class description
464//
465Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
466{
467 Double_t alpha, energy, theta;
468 Double_t size=-1;
469
470 if (fMatrix)
471 {
472 alpha = fMap[0]<0 ? GetVal() : (*fMatrix)[fMap[0]];
473 energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
474 size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
475 //<0 ? 1000 : (*fMatrix)[fMap[1]];
476 theta = 0;
477
478 if (energy<0)
479 energy=size;
480 if (size<0)
481 size=energy;
482
483 if (energy<0 && size<0)
484 energy = size = 1000;
485 }
486 else
487 {
488 alpha = GetVal();
489
490 if (fHillas)
491 size = fHillas->GetSize();
492 energy = fEnergy ? fEnergy->GetVal() : (fHillas?fHillas->GetSize():1000);
493 theta = fPointPos ? fPointPos->GetZd() : 0;
494 }
495
496 // enhance histogram if necessary
497 UpdateAlphaTime();
498
499 // Fill histograms
500 const Int_t ibin = fHist.Fill(theta, energy, TMath::Abs(alpha), w);
501 SetBin(ibin);
502
503 if (!fSkipHistTime)
504 fHistTime.Fill(TMath::Abs(alpha), w);
505
506 return kTRUE;
507}
508
509// --------------------------------------------------------------------------
510//
511// Paint the integral and the error on top of the histogram
512//
513void MHAlpha::PaintText(Double_t val, Double_t error) const
514{
515 TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
516 text.SetBit(TLatex::kTextNDC);
517 text.SetTextSize(0.04);
518 text.Paint();
519}
520
521// --------------------------------------------------------------------------
522//
523// Paint the integral and the error on top of the histogram
524//
525void MHAlpha::PaintText(const TH1D &h) const
526{
527 Double_t sumv = 0;
528 Double_t sume = 0;
529
530 for (int i=0; i<h.GetNbinsX(); i++)
531 {
532 sumv += h.GetBinContent(i+1);
533 sume += h.GetBinError(i+1);
534 }
535 PaintText(sumv, sume);
536}
537// --------------------------------------------------------------------------
538//
539// Update the projections
540//
541void MHAlpha::Paint(Option_t *opt)
542{
543 // Note: Any object cannot be added to a pad twice!
544 // The object is searched by gROOT->FindObject only in
545 // gPad itself!
546 TVirtualPad *padsave = gPad;
547
548 TH1D *h0=0;
549
550 TString o(opt);
551
552 if (o==(TString)"proj")
553 {
554 TPaveStats *st=0;
555 for (int x=0; x<4; x++)
556 {
557 TVirtualPad *p=padsave->GetPad(x+1);
558 if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
559 continue;
560
561 if (st->GetOptStat()==11)
562 continue;
563
564 const Double_t y1 = st->GetY1NDC();
565 const Double_t y2 = st->GetY2NDC();
566 const Double_t x1 = st->GetX1NDC();
567 const Double_t x2 = st->GetX2NDC();
568
569 st->SetY1NDC((y2-y1)/3+y1);
570 st->SetX1NDC((x2-x1)/3+x1);
571 st->SetOptStat(11);
572 }
573
574 padsave->cd(1);
575
576 TH1D *hon = (TH1D*)gPad->FindObject("Proj");
577 if (hon)
578 {
579 TH1D *dum = fHist.ProjectionZ("dum", 0, fHist.GetNbinsX()+1, 0, fHist.GetNbinsY()+1);
580 dum->SetDirectory(0);
581 hon->Reset();
582 hon->Add(dum);
583 delete dum;
584
585 if (fOffData)
586 {
587 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
588 if (hoff)
589 {
590 TH1D *dum = fOffData->ProjectionZ("dum", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1);
591 dum->SetDirectory(0);
592 hoff->Reset();
593 hoff->Add(dum);
594 delete dum;
595
596 const Double_t alpha = fFit.Scale(*hoff, *hon);
597
598 hon->SetMaximum();
599 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
600
601 // BE CARFEULL: This is a really weird workaround!
602 hoff->SetMaximum(alpha);
603
604 // For some reason the line-color is resetted
605 hoff->SetLineColor(kRed);
606
607 if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
608 {
609 h0->Reset();
610 h0->Add(hoff, hon, -1);
611 const Float_t min = h0->GetMinimum()*1.05;
612 hon->SetMinimum(min<0 ? min : 0);
613 }
614 }
615 }
616 else
617 hon->SetMinimum(0);
618 }
619 FitEnergyBins();
620 FitThetaBins();
621 }
622
623 if (o==(TString)"variable")
624 if ((h0 = (TH1D*)gPad->FindObject("Proj")))
625 {
626 // Check whether we are dealing with on-off analysis
627 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
628
629 // BE CARFEULL: This is a really weird workaround!
630 const Double_t scale = !hoff || hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();
631
632 // Do not store the 'final' result in fFit
633 MAlphaFitter fit(fFit);
634 fit.Fit(*h0, hoff, scale, kTRUE);
635 fit.PaintResult();
636 }
637
638 if (o==(TString)"time")
639 PaintText(fHTime);
640
641 if (o==(TString)"theta")
642 {
643 TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName()));
644 if (h)
645 {
646 TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
647 h2->SetDirectory(0);
648 h2->Scale(fHTheta.Integral()/h2->Integral());
649 h->Reset();
650 h->Add(h2);
651 delete h2;
652 }
653 PaintText(fHTheta);
654 }
655
656 if (o==(TString)"energy")
657 {
658 TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName()));
659 if (h)
660 {
661 TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
662 h2->SetDirectory(0);
663 h2->Scale(fHEnergy.Integral()/h2->Integral());
664 h->Reset();
665 h->Add(h2);
666 delete h2;
667 }
668 PaintText(fHEnergy);
669
670 if (fHEnergy.GetMaximum()>1)
671 {
672 gPad->SetLogx();
673 gPad->SetLogy();
674 }
675 }
676
677 gPad = padsave;
678}
679
680// --------------------------------------------------------------------------
681//
682// Draw the histogram
683//
684void MHAlpha::Draw(Option_t *opt)
685{
686 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
687
688 /*
689 if (TString(opt).Contains("sizebins", TString::kIgnoreCase))
690 {
691 AppendPad("sizebins");
692 return;
693 }
694 */
695
696 // Do the projection before painting the histograms into
697 // the individual pads
698 AppendPad("proj");
699
700 pad->SetBorderMode(0);
701 pad->Divide(2,2);
702
703 TH1D *h=0;
704
705 pad->cd(1);
706 gPad->SetBorderMode(0);
707
708 h = fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, 0, fHist.GetNbinsY()+1, "E");
709 h->SetBit(TH1::kNoTitle);
710 h->SetStats(kTRUE);
711 h->SetXTitle(fHist.GetZaxis()->GetTitle());
712 h->SetYTitle("Counts");
713 h->SetDirectory(NULL);
714 h->SetMarkerStyle(0);
715 h->SetBit(kCanDelete);
716 h->Draw("");
717
718 if (fOffData)
719 {
720 // To get green on-data
721 //h->SetMarkerColor(kGreen);
722 //h->SetLineColor(kGreen);
723
724 h = fOffData->ProjectionZ("ProjOff", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E");
725 h->SetBit(TH1::kNoTitle);
726 h->SetXTitle(fHist.GetZaxis()->GetTitle());
727 h->SetYTitle("Counts");
728 h->SetDirectory(NULL);
729 h->SetMarkerStyle(0);
730 h->SetBit(kCanDelete);
731 h->SetMarkerColor(kRed);
732 h->SetLineColor(kRed);
733 //h->SetFillColor(18);
734 h->Draw("same"/*"bar same"*/);
735
736 // This is the only way to make it work...
737 // Clone and copy constructor give strange crashes :(
738 h = fOffData->ProjectionZ("ProjOnOff", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E");
739 h->SetBit(TH1::kNoTitle);
740 h->SetXTitle(fHist.GetZaxis()->GetTitle());
741 h->SetYTitle("Counts");
742 h->SetDirectory(NULL);
743 h->SetMarkerStyle(0);
744 h->SetBit(kCanDelete);
745 h->SetMarkerColor(kBlue);
746 h->SetLineColor(kBlue);
747 h->Draw("same");
748
749 TLine lin;
750 lin.SetLineStyle(kDashed);
751 lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
752 }
753
754 // After the Alpha-projection has been drawn. Fit the histogram
755 // and paint the result into this pad
756 AppendPad("variable");
757
758 if (fHEnergy.GetNbinsX()>1 || fHEnergy.GetBinContent(1)>0)
759 {
760 pad->cd(2);
761 gPad->SetBorderMode(0);
762 gPad->SetGridx();
763 gPad->SetGridy();
764 fHEnergy.Draw();
765
766 AppendPad("energy");
767
768 h = (TH1D*)fHist.Project3D("y");
769 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
770 h->SetXTitle("E [GeV]");
771 h->SetYTitle("Counts");
772 h->SetDirectory(NULL);
773 h->SetMarkerStyle(kFullDotSmall);
774 h->SetBit(kCanDelete);
775 h->SetMarkerColor(kCyan);
776 h->SetLineColor(kCyan);
777 h->Draw("Psame");
778 }
779 else
780 delete pad->GetPad(2);
781
782 if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1 || fHTime.GetBinError(1)>0)
783 {
784 pad->cd(3);
785 gPad->SetBorderMode(0);
786 gPad->SetGridx();
787 gPad->SetGridy();
788 fHTime.Draw();
789 AppendPad("time");
790 }
791 else
792 delete pad->GetPad(3);
793
794 if (fHTheta.GetNbinsX()>1 || fHTheta.GetBinContent(1)>0)
795 {
796 pad->cd(4);
797 gPad->SetGridx();
798 gPad->SetGridy();
799 gPad->SetBorderMode(0);
800 fHTheta.Draw();
801
802 AppendPad("theta");
803
804 h = (TH1D*)fHist.Project3D("x");
805 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
806 h->SetXTitle("\\theta [\\circ]");
807 h->SetYTitle("Counts");
808 h->SetDirectory(NULL);
809 h->SetMarkerStyle(kFullDotSmall);
810 h->SetBit(kCanDelete);
811 h->SetMarkerColor(kCyan);
812 h->SetLineColor(kCyan);
813 h->Draw("Psame");
814 }
815 else
816 delete pad->GetPad(4);
817}
818
819void MHAlpha::DrawAll(Bool_t newc)
820{
821 if (!newc && !fDisplay)
822 return;
823
824 // FIXME: Do in Paint if special option given!
825 TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("SizeBins");
826 Int_t n = fHist.GetNbinsY();
827 Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
828 c.Divide(nc, nc, 1e-10, 1e-10);
829 gPad->SetBorderMode(0);
830 gPad->SetFrameBorderMode(0);
831
832 // Do not store the 'final' result in fFit
833 MAlphaFitter fit(fFit);
834
835 for (int i=1; i<=fHist.GetNbinsY(); i++)
836 {
837 c.cd(i);
838 gPad->SetBorderMode(0);
839 gPad->SetFrameBorderMode(0);
840
841 TH1D *hon = fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, i, i, "E");
842 hon->SetBit(TH1::kNoTitle);
843 hon->SetStats(kTRUE);
844 hon->SetXTitle(fHist.GetZaxis()->GetTitle());
845 hon->SetYTitle("Counts");
846 hon->SetDirectory(NULL);
847 hon->SetMarkerStyle(0);
848 hon->SetBit(kCanDelete);
849 hon->Draw("");
850
851 TH1D *hof = 0;
852 Double_t alpha = 1;
853
854 if (fOffData)
855 {
856 hon->SetMarkerColor(kGreen);
857
858 hof = fOffData->ProjectionZ("ProjOff", 0, fOffData->GetNbinsX()+1, i, i, "E");
859 hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
860 hof->SetXTitle(fHist.GetZaxis()->GetTitle());
861 hof->SetYTitle("Counts");
862 hof->SetDirectory(NULL);
863 hof->SetMarkerStyle(0);
864 hof->SetBit(kCanDelete);
865 hof->SetMarkerColor(kRed);
866 hof->SetLineColor(kRed);
867 hof->Draw("same");
868
869 alpha = fit.Scale(*hof, *hon);
870
871 hon->SetMaximum();
872 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
873
874 TH1D *diff = new TH1D(*hon);
875 diff->Add(hof, -1);
876 diff->SetBit(TH1::kNoTitle);
877 diff->SetXTitle(fHist.GetZaxis()->GetTitle());
878 diff->SetYTitle("Counts");
879 diff->SetDirectory(NULL);
880 diff->SetMarkerStyle(0);
881 diff->SetBit(kCanDelete);
882 diff->SetMarkerColor(kBlue);
883 diff->Draw("same");
884
885 TLine lin;
886 lin.SetLineStyle(kDashed);
887 lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
888
889 const Float_t min = diff->GetMinimum()*1.05;
890 hon->SetMinimum(min<0 ? min : 0);
891 }
892
893 if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
894 {
895 *fLog << dbg << "Bin " << i << ": sigmaexc=" << fit.GetEventsExcess()/fit.GetErrorExcess() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
896 fit.DrawResult();
897 }
898 /*
899 if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
900 {
901 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
902 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
903 }*/
904 }
905
906}
907
908void MHAlpha::DrawNicePlot(Bool_t newc, const char *title, const char *watermark, Int_t binlo, Int_t binhi)
909{
910 if (!newc && !fDisplay)
911 return;
912
913 if (!fOffData)
914 return;
915
916 // Open and setup canvas/pad
917 TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("ThetsSq");
918
919 c.SetBorderMode(0);
920 c.SetFrameBorderMode(0);
921 c.SetFillColor(kWhite);
922
923 c.SetLeftMargin(0.12);
924 c.SetRightMargin(0.01);
925 c.SetBottomMargin(0.16);
926 c.SetTopMargin(0.18);
927
928 c.SetGridy();
929
930 gStyle->SetOptStat(0);
931
932 // Get on-data
933 TH1D *hon = (TH1D*)fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, binlo, binhi, "E");
934 hon->SetDirectory(NULL);
935 hon->SetBit(kCanDelete);
936 hon->SetMarkerSize(0);
937 hon->SetLineWidth(2);
938 hon->SetLineColor(kBlack);
939 hon->SetMarkerColor(kBlack);
940
941 // Get off-data
942 TH1D *hoff = 0;
943 if (fOffData)
944 {
945 hoff = (TH1D*)fOffData->ProjectionZ("ProjOff", 0, fHist.GetNbinsX()+1, binlo, binhi, "E");
946 hoff->SetDirectory(NULL);
947 hoff->SetBit(kCanDelete);
948 hoff->SetFillColor(17);
949 hoff->SetMarkerSize(0);
950 hoff->SetLineColor(kBlack);
951 hoff->SetMarkerColor(kBlack);
952 }
953
954 // Setup plot which is drawn first
955 TH1D *h = hoff ? hoff : hon;
956 h->GetXaxis()->SetLabelSize(0.06);
957 h->GetXaxis()->SetTitleSize(0.06);
958 h->GetXaxis()->SetTitleOffset(0.95);
959 h->GetYaxis()->SetLabelSize(0.06);
960 h->GetYaxis()->SetTitleSize(0.06);
961 h->GetYaxis()->CenterTitle();
962 h->SetYTitle("Counts");
963 h->SetTitleSize(0.07);
964 h->SetTitle("");
965
966 const Double_t imax = fFit.GetSignalIntegralMax();
967 if (imax<1)
968 h->GetXaxis()->SetRangeUser(0, 0.6*0.6);
969
970 // scale off-data
971
972 MAlphaFitter fit(fFit);
973 fit.ScaleAndFit(*hon, hoff);
974
975 hon->SetMinimum(0);
976 hoff->SetMinimum(0);
977
978 // draw data
979 if (hoff)
980 {
981 hoff->SetMaximum(TMath::Max(hon->GetMaximum(),hoff->GetMaximum())*1.1);
982 hoff->Draw("bar");
983 hon->Draw("same");
984 }
985 else
986 {
987 hon->SetMaximum();
988 hon->Draw();
989 }
990
991 // draw a preliminary tag
992 TLatex text;
993 text.SetTextColor(kWhite);
994 text.SetBit(TLatex::kTextNDC);
995 text.SetTextSize(0.07);
996 text.SetTextAngle(2.5);
997
998 TString wm(watermark);
999 if (binlo>=1 || binhi<hon->GetNbinsX())
1000 {
1001 wm += wm.IsNull() ? "(" : " (";
1002 if (binlo>=1)
1003 wm += Form("%.1fGeV", fHist.GetYaxis()->GetBinLowEdge(binlo));
1004 wm += "-";
1005 if (binhi<hon->GetNbinsX())
1006 wm += Form("%.1fGeV", fHist.GetYaxis()->GetBinLowEdge(binhi+1));
1007 wm += ")";
1008 }
1009 if (!wm.IsNull())
1010 text.DrawLatex(0.45, 0.2, wm);
1011 //enum { kTextNDC = BIT(14) };
1012
1013 // draw line showing cut
1014 TLine line;
1015 line.SetLineColor(14);
1016 line.SetLineStyle(7);
1017 line.DrawLine(imax, 0, imax, h->GetMaximum());
1018
1019 // Add a title above the plot
1020 TPaveText *pave=new TPaveText(0.12, 0.83, 0.99, 0.975, "blNDC");
1021 pave->SetBorderSize(1);
1022 pave->SetLabel(title);
1023
1024 char txt[1000];
1025 TText *ptxt;
1026 sprintf(txt, " ");
1027 ptxt = pave->AddText(txt);
1028 ptxt->SetTextAlign(23);
1029
1030 //sprintf(txt, "Significance %.1f\\sigma, off-scale %.2f (\\omega=%.2f\\circ)",
1031 // fFit.GetSignificance(), fFit.GetScaleFactor(), fFit.GetGausSigma());
1032 sprintf(txt, "Significance %.1f\\sigma, off-scale %.2f",
1033 fit.GetSignificance(), fit.GetScaleFactor());
1034 ptxt = pave->AddText(txt);
1035 ptxt->SetTextAlign(23);
1036
1037 sprintf(txt, "%.1f excess events, %.1f background events",
1038 fit.GetEventsExcess(), fit.GetEventsBackground());
1039 ptxt = pave->AddText(txt);
1040 ptxt->SetTextAlign(23);
1041 pave->SetBit(kCanDelete);
1042 pave->Draw();
1043}
1044
1045Bool_t MHAlpha::Finalize()
1046{
1047 if (!FitAlpha())
1048 {
1049 *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
1050 return kTRUE;
1051 }
1052
1053 // Store the final result in fFit
1054 *fLog << all;
1055 fFit.Print("result");
1056
1057 fResult->SetVal(fFit.GetMinimizationValue());
1058 fSigma->SetVal(fFit.GetGausSigma());
1059
1060 if (!fSkipHistEnergy)
1061 {
1062 *fLog << inf << "Processing energy bins..." << endl;
1063 FitEnergyBins();
1064 }
1065 if (!fSkipHistTheta)
1066 {
1067 *fLog << inf << "Processing theta bins..." << endl;
1068 FitThetaBins();
1069 }
1070 if (!fSkipHistTime)
1071 {
1072 *fLog << inf << "Processing time bins..." << endl;
1073 UpdateAlphaTime(kTRUE);
1074 MH::RemoveFirstBin(fHTime);
1075 }
1076
1077 if (fOffData)
1078 DrawAll(kFALSE);
1079
1080 return kTRUE;
1081}
1082
1083// --------------------------------------------------------------------------
1084//
1085// You can use this function if you want to use a MHMatrix instead of
1086// MMcEvt. This function adds all necessary columns to the
1087// given matrix. Afterward you should fill the matrix with the corresponding
1088// data (eg from a file by using MHMatrix::Fill). If you now loop
1089// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
1090// will take the values from the matrix instead of the containers.
1091//
1092// It takes fSkipHist* into account!
1093//
1094void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
1095{
1096 if (fMatrix)
1097 return;
1098
1099 fMatrix = mat;
1100
1101 fMap[0] = fMatrix->AddColumn(GetParameterRule());
1102 fMap[1] = -1;
1103 fMap[2] = -1;
1104 fMap[3] = -1;
1105 fMap[4] = -1;
1106
1107 if (!fSkipHistEnergy)
1108 if (type==0)
1109 {
1110 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
1111 fMap[2] = -1;
1112 }
1113 else
1114 {
1115 fMap[1] = -1;
1116 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
1117 }
1118
1119 if (!fSkipHistTheta)
1120 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
1121
1122 // if (!fSkipHistTime)
1123 // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
1124}
1125
1126void MHAlpha::StopMapping()
1127{
1128 fMatrix = NULL;
1129}
1130
1131void MHAlpha::ApplyScaling()
1132{
1133 if (!fOffData)
1134 return;
1135
1136 fFit.ApplyScaling(fHist, *const_cast<TH3D*>(fOffData));
1137}
1138
1139Int_t MHAlpha::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1140{
1141 Bool_t rc = kFALSE;
1142 if (IsEnvDefined(env, prefix, "NumTimeBins", print))
1143 {
1144 SetNumTimeBins(GetEnvValue(env, prefix, "NumTimeBins", fNumTimeBins));
1145 rc = kTRUE;
1146 }
1147 if (IsEnvDefined(env, prefix, "ForceUsingSize", print))
1148 {
1149 fForceUsingSize = GetEnvValue(env, prefix, "ForceUsingSize", fForceUsingSize);
1150 rc = kTRUE;
1151 }
1152 return rc;
1153}
1154
1155Int_t MHAlpha::DistancetoPrimitive(Int_t px, Int_t py)
1156{
1157 // If pad has no subpad return (we are in one of the subpads)
1158 if (gPad->GetPad(1)==NULL)
1159 return 9999;
1160
1161 // If pad has a subpad we are in the main pad. Check its value.
1162 return gPad->GetPad(1)->DistancetoPrimitive(px,py)==0 ? 0 : 9999;
1163}
Note: See TracBrowser for help on using the repository browser.