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

Last change on this file since 7764 was 7755, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 26.3 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 "MGeomCam.h"
57#include "MSrcPosCam.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);
74ClassImp(MAlphaFitter);
75
76using namespace std;
77
78// --------------------------------------------------------------------------
79//
80// Default Constructor
81//
82MHAlpha::MHAlpha(const char *name, const char *title)
83 : fNameParameter("MHillasSrc"), fParameter(0),
84 fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0),
85 fPointPos(0), fTimeEffOn(0), fTime(0), fNumTimeBins(10),
86 fHillas(0), fMatrix(0), fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE),
87 fSkipHistEnergy(kFALSE), fForceUsingSize(kFALSE)
88{
89 //
90 // set the name and title of this object
91 //
92 fName = name ? name : "MHAlpha";
93 fTitle = title ? title : "Alpha plot";
94
95 fHist.SetName("Alpha");
96 fHist.SetTitle("Alpha");
97 fHist.SetXTitle("\\Theta [deg]");
98 //fHist.SetYTitle("E_{est} [GeV]");
99 fHist.SetZTitle("|\\alpha| [\\circ]");
100 fHist.SetDirectory(NULL);
101 fHist.UseCurrentStyle();
102
103 // Main histogram
104 fHistTime.SetName("Alpha");
105 fHistTime.SetXTitle("|\\alpha| [\\circ]");
106 fHistTime.SetYTitle("Counts");
107 fHistTime.UseCurrentStyle();
108 fHistTime.SetDirectory(NULL);
109
110
111 fHEnergy.SetName("Excess");
112 //fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
113 //fHEnergy.SetXTitle("E_{est} [GeV]");
114 fHEnergy.SetYTitle("N_{exc}");
115 fHEnergy.SetDirectory(NULL);
116 fHEnergy.UseCurrentStyle();
117
118 fHTheta.SetName("ExcessTheta");
119 fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
120 fHTheta.SetXTitle("\\Theta [\\circ]");
121 fHTheta.SetYTitle("N_{exc}");
122 fHTheta.SetDirectory(NULL);
123 fHTheta.UseCurrentStyle();
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.Sumw2();
137
138 MBinning binsa, binse, binst;
139 binsa.SetEdges(18, 0, 90);
140 binse.SetEdgesLog(15, 10, 100000);
141 binst.SetEdgesASin(67, -0.005, 0.665);
142 binse.Apply(fHEnergy);
143 binst.Apply(fHTheta);
144 binsa.Apply(fHistTime);
145
146 MH::SetBinning(&fHist, &binst, &binse, &binsa);
147}
148
149Float_t MHAlpha::FitEnergyBins(Bool_t paint)
150{
151 // Do not store the 'final' result in fFit
152 MAlphaFitter fit(fFit);
153
154 const Int_t n = fHist.GetNbinsY();
155
156 fHEnergy.SetEntries(0);
157
158 Float_t mean = 0;
159 Int_t num = 0;
160 for (int i=1; i<=n; i++)
161 {
162 if (fit.FitEnergy(fHist, fOffData, i))
163 {
164 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
165 if (fit.GetSignificanceExc()>1)
166 fHEnergy.SetBinError(i, fit.GetEventsExcess()/fit.GetSignificanceExc());
167 else
168 fHEnergy.SetBinError(i, fit.GetEventsExcess());
169
170 if (fit.GetSignificanceExc()>1)
171 {
172 mean += fit.GetSignificanceExc()*fit.GetSignificanceExc();
173 num++;
174 }
175 }
176 }
177 return TMath::Sqrt(mean)/num;
178}
179
180void MHAlpha::FitThetaBins(Bool_t paint)
181{
182 // Do not store the 'final' result in fFit
183 MAlphaFitter fit(fFit);
184
185 const Int_t n = fHist.GetNbinsX();
186
187 fHTheta.SetEntries(0);
188
189 for (int i=1; i<=n; i++)
190 {
191 if (fit.FitTheta(fHist, fOffData, i))
192 {
193 fHTheta.SetBinContent(i, fit.GetEventsExcess());
194 if (fit.GetSignificanceExc()>1)
195 fHTheta.SetBinError(i, fit.GetEventsExcess()/fit.GetSignificanceExc());
196 else
197 fHTheta.SetBinError(i, fit.GetEventsExcess());
198 }
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
298 //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
299 //if (!fExcess)
300 // return kFALSE;
301
302 fLastTime = MTime();
303 fNumRebin = fNumTimeBins;
304
305 MBinning binst, binse, binsa;
306 binst.SetEdges(fOffData ? *fOffData : fHist, 'x');
307 binse.SetEdges(fOffData ? *fOffData : fHist, 'y');
308 binsa.SetEdges(fOffData ? *fOffData : fHist, 'z');
309 if (!fOffData)
310 {
311 if (!fPointPos)
312 binst.SetEdges(1, 0, 60);
313 else
314 binst.SetEdges(*pl, "BinningTheta");
315
316 if (!fEnergy && !fHillas)
317 binse.SetEdges(1, 10, 100000);
318 else
319 if (fEnergy)
320 binse.SetEdges(*pl, "BinningEnergyEst");
321 else
322 binse.SetEdges(*pl, "BinningSize");
323
324 binsa.SetEdges(*pl, Form("Binning%s", ClassName()+2));
325 }
326 else
327 {
328 fHEnergy.SetTitle(fOffData->GetTitle());
329 fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle());
330 fHist.SetYTitle(fOffData->GetYaxis()->GetTitle());
331 }
332
333 binse.Apply(fHEnergy);
334 binst.Apply(fHTheta);
335 binsa.Apply(fHistTime);
336 MH::SetBinning(&fHist, &binst, &binse, &binsa);
337
338 MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
339 if (!fit)
340 *fLog << warn << "MAlphaFitter not found... ignored." << endl;
341 else
342 fit->Copy(fFit);
343
344 *fLog << inf;
345 fFit.Print();
346
347 return kTRUE;
348}
349
350void MHAlpha::InitAlphaTime(const MTime &t)
351{
352 //
353 // If this is the first call we have to initialize the time-histogram
354 //
355 MBinning bins;
356 bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
357 bins.Apply(fHTime);
358
359 fLastTime=t;
360}
361
362void MHAlpha::UpdateAlphaTime(Bool_t final)
363{
364 if (!fTimeEffOn)
365 return;
366
367 if (!final)
368 {
369 if (fLastTime==MTime() && *fTimeEffOn!=MTime())
370 {
371 InitAlphaTime(*fTimeEffOn);
372 return;
373 }
374
375 if (fLastTime!=*fTimeEffOn)
376 {
377 fLastTime=*fTimeEffOn;
378 fNumRebin--;
379 }
380
381 if (fNumRebin>0)
382 return;
383 }
384
385 // Check new 'last time'
386 MTime *time = final ? fTime : fTimeEffOn;
387
388 if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
389 {
390 *fLog << warn << "WARNING - New time-stamp " << *time << " lower" << endl;
391 *fLog << "than upper edge of histogram... skipped." << endl;
392 *fLog << "This should not happen. Maybe you started you eventloop with" << endl;
393 *fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
394 fNumRebin++;
395 return;
396 }
397
398 // Fit time histogram
399 MAlphaFitter fit(fFit);
400
401 TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
402 const Bool_t rc = fit.ScaleAndFit(fHistTime, h);
403 if (h)
404 delete h;
405
406 if (!rc)
407 return;
408
409 // Reset Histogram
410 fHistTime.Reset();
411 fHistTime.SetEntries(0);
412
413 //
414 // Prepare Histogram
415 //
416 if (final)
417 time->Plus1ns();
418
419 // Get number of bins
420 const Int_t n = fHTime.GetNbinsX();
421
422 // Enhance binning
423 MBinning bins;
424 bins.SetEdges(fHTime, 'x');
425 bins.AddEdge(time->GetAxisTime());
426 bins.Apply(fHTime);
427
428 //
429 // Fill histogram
430 //
431 fHTime.SetBinContent(n+1, fit.GetEventsExcess());
432 if (fit.GetSignificanceExc()>1)
433 fHTime.SetBinError(n+1, fit.GetEventsExcess()/fit.GetSignificanceExc());
434 else
435 fHTime.SetBinError(n+1, fit.GetEventsExcess());
436
437 *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl;
438
439 fNumRebin = fNumTimeBins;
440}
441
442// --------------------------------------------------------------------------
443//
444// Fill the histogram. For details see the code or the class description
445//
446Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
447{
448 Double_t alpha, energy, theta;
449 Double_t size=-1;
450
451 if (fMatrix)
452 {
453 alpha = fMap[0]<0 ? GetVal() : (*fMatrix)[fMap[0]];
454 energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
455 size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
456 //<0 ? 1000 : (*fMatrix)[fMap[1]];
457 theta = 0;
458
459 if (energy<0)
460 energy=size;
461 if (size<0)
462 size=energy;
463
464 if (energy<0 && size<0)
465 energy = size = 1000;
466 }
467 else
468 {
469 alpha = GetVal();
470
471 if (fHillas)
472 size = fHillas->GetSize();
473 energy = fEnergy ? fEnergy->GetVal() : (fHillas?fHillas->GetSize():1000);
474 theta = fPointPos ? fPointPos->GetZd() : 0;
475 }
476
477 // enhance histogram if necessary
478 UpdateAlphaTime();
479
480 // Fill histograms
481 fHist.Fill(theta, energy, TMath::Abs(alpha), w);
482
483 if (!fSkipHistTime)
484 fHistTime.Fill(TMath::Abs(alpha), w);
485
486 return kTRUE;
487}
488
489// --------------------------------------------------------------------------
490//
491// Paint the integral and the error on top of the histogram
492//
493void MHAlpha::PaintText(Double_t val, Double_t error) const
494{
495 TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
496 text.SetBit(TLatex::kTextNDC);
497 text.SetTextSize(0.04);
498 text.Paint();
499}
500
501// --------------------------------------------------------------------------
502//
503// Paint the integral and the error on top of the histogram
504//
505void MHAlpha::PaintText(const TH1D &h) const
506{
507 Double_t sumv = 0;
508 Double_t sume = 0;
509
510 for (int i=0; i<h.GetNbinsX(); i++)
511 {
512 sumv += h.GetBinContent(i+1);
513 sume += h.GetBinError(i+1);
514 }
515 PaintText(sumv, sume);
516}
517// --------------------------------------------------------------------------
518//
519// Update the projections
520//
521void MHAlpha::Paint(Option_t *opt)
522{
523 // Note: Any object cannot be added to a pad twice!
524 // The object is searched by gROOT->FindObject only in
525 // gPad itself!
526 TVirtualPad *padsave = gPad;
527
528 TH1D *h0=0;
529
530 TString o(opt);
531 if (o==(TString)"proj")
532 {
533 TPaveStats *st=0;
534 for (int x=0; x<4; x++)
535 {
536 TVirtualPad *p=padsave->GetPad(x+1);
537 if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
538 continue;
539
540 if (st->GetOptStat()==11)
541 continue;
542
543 const Double_t y1 = st->GetY1NDC();
544 const Double_t y2 = st->GetY2NDC();
545 const Double_t x1 = st->GetX1NDC();
546 const Double_t x2 = st->GetX2NDC();
547
548 st->SetY1NDC((y2-y1)/3+y1);
549 st->SetX1NDC((x2-x1)/3+x1);
550 st->SetOptStat(11);
551 }
552
553 padsave->cd(1);
554 TH1D *hon = fHist.ProjectionZ("Proj");
555 if (fOffData)
556 {
557 TH1D *hoff = fOffData->ProjectionZ("ProjOff");
558 const Double_t alpha = fFit.Scale(*hoff, *hon);
559
560 hon->SetMaximum();
561 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
562
563 // BE CARFEULL: This is a really weird workaround!
564 hoff->SetMaximum(alpha);
565
566 if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
567 {
568 h0->Reset();
569 h0->Add(hoff, hon, -1);
570 const Float_t min = h0->GetMinimum()*1.05;
571 hon->SetMinimum(min<0 ? min : 0);
572 }
573 }
574 else
575 hon->SetMinimum(0);
576 FitEnergyBins();
577 FitThetaBins();
578 }
579
580 if (o==(TString)"variable")
581 if ((h0 = (TH1D*)gPad->FindObject("Proj")))
582 {
583 // Check whether we are dealing with on-off analysis
584 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
585
586 // BE CARFEULL: This is a really weird workaround!
587 const Double_t scale = !hoff || hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();
588
589 // Do not store the 'final' result in fFit
590 MAlphaFitter fit(fFit);
591 fit.Fit(*h0, hoff, scale, kTRUE);
592 fit.PaintResult();
593 }
594
595 if (o==(TString)"time")
596 PaintText(fHTime);//.Integral(), 0);
597
598 if (o==(TString)"theta")
599 {
600 TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName()));
601 if (h)
602 {
603 TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
604 h2->SetDirectory(0);
605 h2->Scale(fHTheta.Integral()/h2->Integral());
606 h->Reset();
607 h->Add(h2);
608 delete h2;
609 }
610 PaintText(fHTheta);//.Integral(), 0);
611 }
612
613 if (o==(TString)"energy")
614 {
615 TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName()));
616 if (h)
617 {
618 TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
619 h2->SetDirectory(0);
620 h2->Scale(fHEnergy.Integral()/h2->Integral());
621 h->Reset();
622 h->Add(h2);
623 delete h2;
624 }
625 PaintText(fHEnergy);//.Integral(), 0);
626
627 if (fHEnergy.GetMaximum()>1)
628 {
629 gPad->SetLogx();
630 gPad->SetLogy();
631 }
632 }
633
634 gPad = padsave;
635}
636
637// --------------------------------------------------------------------------
638//
639// Draw the histogram
640//
641void MHAlpha::Draw(Option_t *opt)
642{
643 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
644
645 /*
646 if (TString(opt).Contains("sizebins", TString::kIgnoreCase))
647 {
648 AppendPad("sizebins");
649 return;
650 }
651 */
652
653 // Do the projection before painting the histograms into
654 // the individual pads
655 AppendPad("proj");
656
657 pad->SetBorderMode(0);
658 pad->Divide(2,2);
659
660 TH1D *h=0;
661
662 pad->cd(1);
663 gPad->SetBorderMode(0);
664
665 h = fHist.ProjectionZ("Proj", -1, 9999, -1, 9999, "E");
666 h->SetBit(TH1::kNoTitle);
667 h->SetStats(kTRUE);
668 h->SetXTitle(fHist.GetZaxis()->GetTitle());
669 h->SetYTitle("Counts");
670 h->SetDirectory(NULL);
671 h->SetMarkerStyle(kFullDotMedium);
672 h->SetBit(kCanDelete);
673 h->Draw("");
674
675 if (fOffData)
676 {
677 h->SetMarkerColor(kGreen);
678
679 h = fOffData->ProjectionZ("ProjOff", -1, 9999, -1, 9999, "E");
680 h->SetBit(TH1::kNoTitle);
681 h->SetXTitle(fHist.GetZaxis()->GetTitle());
682 h->SetYTitle("Counts");
683 h->SetDirectory(NULL);
684 h->SetMarkerStyle(kFullDotMedium);
685 h->SetBit(kCanDelete);
686 h->SetMarkerColor(kRed);
687 h->Draw("same");
688
689 h = (TH1D*)h->Clone("ProjOnOff");
690 h->SetBit(TH1::kNoTitle);
691 h->SetXTitle(fHist.GetZaxis()->GetTitle());
692 h->SetYTitle("Counts");
693 h->SetDirectory(NULL);
694 h->SetMarkerStyle(kFullDotMedium);
695 h->SetBit(kCanDelete);
696 h->SetMarkerColor(kBlue);
697 h->Draw("same");
698
699 TLine lin;
700 lin.SetLineStyle(kDashed);
701 lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
702 }
703
704 // After the Alpha-projection has been drawn. Fit the histogram
705 // and paint the result into this pad
706 AppendPad("variable");
707
708 if (fHEnergy.GetNbinsX()>1)
709 {
710 pad->cd(2);
711 gPad->SetBorderMode(0);
712 gPad->SetGridx();
713 gPad->SetGridy();
714 fHEnergy.Draw();
715
716 AppendPad("energy");
717
718 h = (TH1D*)fHist.Project3D("y");
719 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
720 h->SetXTitle("E [GeV]");
721 h->SetYTitle("Counts");
722 h->SetDirectory(NULL);
723 h->SetMarkerStyle(kFullDotMedium);
724 h->SetBit(kCanDelete);
725 h->SetMarkerColor(kCyan);
726 h->SetLineColor(kCyan);
727 h->Draw("Psame");
728 }
729 else
730 delete pad->GetPad(2);
731
732 if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
733 {
734 pad->cd(3);
735 gPad->SetBorderMode(0);
736 gPad->SetGridx();
737 gPad->SetGridy();
738 fHTime.Draw();
739 AppendPad("time");
740 }
741 else
742 delete pad->GetPad(3);
743
744 if (fHTheta.GetNbinsX()>1)
745 {
746 pad->cd(4);
747 gPad->SetGridx();
748 gPad->SetGridy();
749 gPad->SetBorderMode(0);
750 fHTheta.Draw();
751
752 AppendPad("theta");
753
754 h = (TH1D*)fHist.Project3D("x");
755 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
756 h->SetXTitle("\\theta [\\circ]");
757 h->SetYTitle("Counts");
758 h->SetDirectory(NULL);
759 h->SetMarkerStyle(kFullDotMedium);
760 h->SetBit(kCanDelete);
761 h->SetMarkerColor(kCyan);
762 h->SetLineColor(kCyan);
763 h->Draw("Psame");
764 }
765 else
766 delete pad->GetPad(4);
767}
768
769void MHAlpha::DrawAll(Bool_t newc)
770{
771 if (newc && !fDisplay)
772 return;
773
774 // FIXME: Do in Paint if special option given!
775 TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("SizeBins");
776 Int_t n = fHist.GetNbinsY();
777 Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
778 c.Divide(nc, nc, 1e-10, 1e-10);
779 gPad->SetBorderMode(0);
780 gPad->SetFrameBorderMode(0);
781
782 // Do not store the 'final' result in fFit
783 MAlphaFitter fit(fFit);
784
785 for (int i=1; i<=fHist.GetNbinsY(); i++)
786 {
787 c.cd(i);
788 gPad->SetBorderMode(0);
789 gPad->SetFrameBorderMode(0);
790
791 TH1D *hon = fHist.ProjectionZ("Proj", -1, 9999, i, i, "E");
792 hon->SetBit(TH1::kNoTitle);
793 hon->SetStats(kTRUE);
794 hon->SetXTitle(fHist.GetZaxis()->GetTitle());
795 hon->SetYTitle("Counts");
796 hon->SetDirectory(NULL);
797 hon->SetMarkerStyle(kFullDotMedium);
798 hon->SetBit(kCanDelete);
799 hon->Draw("");
800
801 TH1D *hof = 0;
802 Double_t alpha = 1;
803
804 if (fOffData)
805 {
806 hon->SetMarkerColor(kGreen);
807
808 hof = fOffData->ProjectionZ("ProjOff", -1, 9999, i, i, "E");
809 hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
810 hof->SetXTitle(fHist.GetZaxis()->GetTitle());
811 hof->SetYTitle("Counts");
812 hof->SetDirectory(NULL);
813 hof->SetMarkerStyle(kFullDotMedium);
814 hof->SetBit(kCanDelete);
815 hof->SetMarkerColor(kRed);
816 hof->Draw("same");
817
818 alpha = fit.Scale(*hof, *hon);
819
820 hon->SetMaximum();
821 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
822
823 TH1D *diff = new TH1D(*hon);
824 diff->Add(hof, -1);
825 diff->SetBit(TH1::kNoTitle);
826 diff->SetXTitle(fHist.GetZaxis()->GetTitle());
827 diff->SetYTitle("Counts");
828 diff->SetDirectory(NULL);
829 diff->SetMarkerStyle(kFullDotMedium);
830 diff->SetBit(kCanDelete);
831 diff->SetMarkerColor(kBlue);
832 diff->Draw("same");
833
834 TLine lin;
835 lin.SetLineStyle(kDashed);
836 lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
837
838 const Float_t min = diff->GetMinimum()*1.05;
839 hon->SetMinimum(min<0 ? min : 0);
840 }
841
842 if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
843 {
844 *fLog << dbg << "Bin " << i << ": sigmaexc=" << fit.GetSignificanceExc() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
845 fit.PaintResult();
846 }
847 /*
848 if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
849 {
850 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
851 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
852 }*/
853 }
854
855}
856
857
858Bool_t MHAlpha::Finalize()
859{
860 //TH1D *h = fHist.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
861 //h->SetDirectory(0);
862 //Bool_t rc = fFit.Fit(*h);
863 //delete h;
864
865 if (!fFit.FitAlpha(fHist, fOffData))
866 {
867 *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
868 return kTRUE;
869 }
870
871 // Store the final result in fFit
872 *fLog << all;
873 fFit.Print("result");
874
875 fResult->SetVal(fFit.GetMinimizationValue());
876 fSigma->SetVal(fFit.GetGausSigma());
877
878 if (!fSkipHistEnergy)
879 {
880 *fLog << inf << "Processing energy bins..." << endl;
881 FitEnergyBins();
882 }
883 if (!fSkipHistTheta)
884 {
885 *fLog << inf << "Processing theta bins..." << endl;
886 FitThetaBins();
887 }
888 if (!fSkipHistTime)
889 {
890 *fLog << inf << "Processing time bins..." << endl;
891 UpdateAlphaTime(kTRUE);
892 MH::RemoveFirstBin(fHTime);
893 }
894
895 if (fOffData)
896 DrawAll(kFALSE);
897
898 return kTRUE;
899}
900
901// --------------------------------------------------------------------------
902//
903// You can use this function if you want to use a MHMatrix instead of
904// MMcEvt. This function adds all necessary columns to the
905// given matrix. Afterward you should fill the matrix with the corresponding
906// data (eg from a file by using MHMatrix::Fill). If you now loop
907// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
908// will take the values from the matrix instead of the containers.
909//
910// It takes fSkipHist* into account!
911//
912void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
913{
914 if (fMatrix)
915 return;
916
917 fMatrix = mat;
918
919 fMap[0] = fMatrix->AddColumn(GetParameterRule());
920 fMap[1] = -1;
921 fMap[2] = -1;
922 fMap[3] = -1;
923 fMap[4] = -1;
924
925 if (!fSkipHistEnergy)
926 if (type==0)
927 {
928 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
929 fMap[2] = -1;
930 }
931 else
932 {
933 fMap[1] = -1;
934 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
935 }
936
937 if (!fSkipHistTheta)
938 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
939
940 // if (!fSkipHistTime)
941 // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
942}
943
944void MHAlpha::StopMapping()
945{
946 fMatrix = NULL;
947}
948
949Int_t MHAlpha::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
950{
951 Bool_t rc = kFALSE;
952 if (IsEnvDefined(env, prefix, "NumTimeBins", print))
953 {
954 SetNumTimeBins(GetEnvValue(env, prefix, "NumTimeBins", fNumTimeBins));
955 rc = kTRUE;
956 }
957 if (IsEnvDefined(env, prefix, "ForceUsingSize", print))
958 {
959 fForceUsingSize = GetEnvValue(env, prefix, "ForceUsingSize", fForceUsingSize);
960 rc = kTRUE;
961 }
962 return rc;
963}
964
965Int_t MHAlpha::DistancetoPrimitive(Int_t px, Int_t py)
966{
967 // If pad has no subpad return (we are in one of the subpads)
968 if (gPad->GetPad(1)==NULL)
969 return 9999;
970
971 // If pad has a subpad we are in the main pad. Check its value.
972 return gPad->GetPad(1)->DistancetoPrimitive(px,py)==0 ? 0 : 9999;
973}
Note: See TracBrowser for help on using the repository browser.