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

Last change on this file since 7562 was 7386, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 25.5 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), fHillas(0),
85 fPointPos(0), fTimeEffOn(0), fTime(0), fNumTimeBins(10),
86 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// Update the projections
504//
505void MHAlpha::Paint(Option_t *opt)
506{
507 // Note: Any object cannot be added to a pad twice!
508 // The object is searched by gROOT->FindObject only in
509 // gPad itself!
510 TVirtualPad *padsave = gPad;
511
512 TH1D *h0=0;
513
514 TString o(opt);
515 if (o==(TString)"proj")
516 {
517 TPaveStats *st=0;
518 for (int x=0; x<4; x++)
519 {
520 TVirtualPad *p=padsave->GetPad(x+1);
521 if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
522 continue;
523
524 if (st->GetOptStat()==11)
525 continue;
526
527 const Double_t y1 = st->GetY1NDC();
528 const Double_t y2 = st->GetY2NDC();
529 const Double_t x1 = st->GetX1NDC();
530 const Double_t x2 = st->GetX2NDC();
531
532 st->SetY1NDC((y2-y1)/3+y1);
533 st->SetX1NDC((x2-x1)/3+x1);
534 st->SetOptStat(11);
535 }
536
537 padsave->cd(1);
538 TH1D *hon = fHist.ProjectionZ("Proj");
539 if (fOffData)
540 {
541 TH1D *hoff = fOffData->ProjectionZ("ProjOff");
542 const Double_t alpha = fFit.Scale(*hoff, *hon);
543
544 hon->SetMaximum();
545 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
546
547 // BE CARFEULL: This is a really weird workaround!
548 hoff->SetMaximum(alpha);
549
550 if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
551 {
552 h0->Reset();
553 h0->Add(hoff, hon, -1);
554 const Float_t min = h0->GetMinimum()*1.05;
555 hon->SetMinimum(min<0 ? min : 0);
556 }
557 }
558 else
559 hon->SetMinimum(0);
560 FitEnergyBins();
561 FitThetaBins();
562 }
563
564 if (o==(TString)"variable")
565 if ((h0 = (TH1D*)gPad->FindObject("Proj")))
566 {
567 // Check whether we are dealing with on-off analysis
568 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
569
570 // BE CARFEULL: This is a really weird workaround!
571 const Double_t scale = !hoff || hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();
572
573 // Do not store the 'final' result in fFit
574 MAlphaFitter fit(fFit);
575 fit.Fit(*h0, hoff, scale, kTRUE);
576 fit.PaintResult();
577 }
578
579 if (o==(TString)"time")
580 PaintText(fHTime.Integral(), 0);
581
582 if (o==(TString)"theta")
583 {
584 TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName()));
585 if (h)
586 {
587 TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
588 h2->SetDirectory(0);
589 h2->Scale(fHTheta.Integral()/h2->Integral());
590 h->Reset();
591 h->Add(h2);
592 delete h2;
593 }
594 PaintText(fHTheta.Integral(), 0);
595 }
596
597 if (o==(TString)"energy")
598 {
599 TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName()));
600 if (h)
601 {
602 TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
603 h2->SetDirectory(0);
604 h2->Scale(fHEnergy.Integral()/h2->Integral());
605 h->Reset();
606 h->Add(h2);
607 delete h2;
608 }
609
610 if (fHEnergy.GetMaximum()>1)
611 {
612 gPad->SetLogx();
613 gPad->SetLogy();
614 }
615 }
616
617 gPad = padsave;
618}
619
620// --------------------------------------------------------------------------
621//
622// Draw the histogram
623//
624void MHAlpha::Draw(Option_t *opt)
625{
626 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
627
628 // Do the projection before painting the histograms into
629 // the individual pads
630 AppendPad("proj");
631
632 pad->SetBorderMode(0);
633 pad->Divide(2,2);
634
635 TH1D *h=0;
636
637 pad->cd(1);
638 gPad->SetBorderMode(0);
639
640 h = fHist.ProjectionZ("Proj", -1, 9999, -1, 9999, "E");
641 h->SetBit(TH1::kNoTitle);
642 h->SetStats(kTRUE);
643 h->SetXTitle(fHist.GetZaxis()->GetTitle());
644 h->SetYTitle("Counts");
645 h->SetDirectory(NULL);
646 h->SetMarkerStyle(kFullDotMedium);
647 h->SetBit(kCanDelete);
648 h->Draw("");
649
650 if (fOffData)
651 {
652 h->SetMarkerColor(kGreen);
653
654 h = fOffData->ProjectionZ("ProjOff", -1, 9999, -1, 9999, "E");
655 h->SetBit(TH1::kNoTitle);
656 h->SetXTitle(fHist.GetZaxis()->GetTitle());
657 h->SetYTitle("Counts");
658 h->SetDirectory(NULL);
659 h->SetMarkerStyle(kFullDotMedium);
660 h->SetBit(kCanDelete);
661 h->SetMarkerColor(kRed);
662 h->Draw("same");
663
664 h = (TH1D*)h->Clone("ProjOnOff");
665 h->SetBit(TH1::kNoTitle);
666 h->SetXTitle(fHist.GetZaxis()->GetTitle());
667 h->SetYTitle("Counts");
668 h->SetDirectory(NULL);
669 h->SetMarkerStyle(kFullDotMedium);
670 h->SetBit(kCanDelete);
671 h->SetMarkerColor(kBlue);
672 h->Draw("same");
673
674 TLine lin;
675 lin.SetLineStyle(kDashed);
676 lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
677 }
678
679 // After the Alpha-projection has been drawn. Fit the histogram
680 // and paint the result into this pad
681 AppendPad("variable");
682
683 if (fHEnergy.GetNbinsX()>1)
684 {
685 pad->cd(2);
686 gPad->SetBorderMode(0);
687 gPad->SetGridx();
688 gPad->SetGridy();
689 fHEnergy.Draw();
690
691 AppendPad("energy");
692
693 h = (TH1D*)fHist.Project3D("y");
694 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
695 h->SetXTitle("E [GeV]");
696 h->SetYTitle("Counts");
697 h->SetDirectory(NULL);
698 h->SetMarkerStyle(kFullDotMedium);
699 h->SetBit(kCanDelete);
700 h->SetMarkerColor(kCyan);
701 h->SetLineColor(kCyan);
702 h->Draw("Psame");
703 }
704 else
705 delete pad->GetPad(2);
706
707 if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
708 {
709 pad->cd(3);
710 gPad->SetBorderMode(0);
711 gPad->SetGridx();
712 gPad->SetGridy();
713 fHTime.Draw();
714 AppendPad("time");
715 }
716 else
717 delete pad->GetPad(3);
718
719 if (fHTheta.GetNbinsX()>1)
720 {
721 pad->cd(4);
722 gPad->SetGridx();
723 gPad->SetGridy();
724 gPad->SetBorderMode(0);
725 fHTheta.Draw();
726
727 AppendPad("theta");
728
729 h = (TH1D*)fHist.Project3D("x");
730 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
731 h->SetXTitle("\\theta [\\circ]");
732 h->SetYTitle("Counts");
733 h->SetDirectory(NULL);
734 h->SetMarkerStyle(kFullDotMedium);
735 h->SetBit(kCanDelete);
736 h->SetMarkerColor(kCyan);
737 h->SetLineColor(kCyan);
738 h->Draw("Psame");
739 }
740 else
741 delete pad->GetPad(4);
742}
743
744void MHAlpha::DrawAll()
745{
746 // FIXME: Do in Paint if special option given!
747 TCanvas *c = new TCanvas;
748 Int_t n = fHist.GetNbinsY();
749 Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
750 c->Divide(nc, nc, 1e-10, 1e-10);
751
752 // Do not store the 'final' result in fFit
753 MAlphaFitter fit(fFit);
754
755 for (int i=1; i<=fHist.GetNbinsY(); i++)
756 {
757 c->cd(i);
758
759 TH1D *hon = fHist.ProjectionZ("Proj", -1, 9999, i, i, "E");
760 hon->SetBit(TH1::kNoTitle);
761 hon->SetStats(kTRUE);
762 hon->SetXTitle(fHist.GetZaxis()->GetTitle());
763 hon->SetYTitle("Counts");
764 hon->SetDirectory(NULL);
765 hon->SetMarkerStyle(kFullDotMedium);
766 hon->SetBit(kCanDelete);
767 hon->Draw("");
768
769 TH1D *hof = 0;
770 Double_t alpha = 1;
771
772 if (fOffData)
773 {
774 hon->SetMarkerColor(kGreen);
775
776 hof = fOffData->ProjectionZ("ProjOff", -1, 9999, i, i, "E");
777 hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
778 hof->SetXTitle(fHist.GetZaxis()->GetTitle());
779 hof->SetYTitle("Counts");
780 hof->SetDirectory(NULL);
781 hof->SetMarkerStyle(kFullDotMedium);
782 hof->SetBit(kCanDelete);
783 hof->SetMarkerColor(kRed);
784 hof->Draw("same");
785
786 alpha = fit.Scale(*hof, *hon);
787
788 hon->SetMaximum();
789 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
790
791 TH1D *diff = new TH1D(*hon);
792 diff->Add(hof, -1);
793 diff->SetBit(TH1::kNoTitle);
794 diff->SetXTitle(fHist.GetZaxis()->GetTitle());
795 diff->SetYTitle("Counts");
796 diff->SetDirectory(NULL);
797 diff->SetMarkerStyle(kFullDotMedium);
798 diff->SetBit(kCanDelete);
799 diff->SetMarkerColor(kBlue);
800 diff->Draw("same");
801
802 TLine lin;
803 lin.SetLineStyle(kDashed);
804 lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
805
806 const Float_t min = diff->GetMinimum()*1.05;
807 hon->SetMinimum(min<0 ? min : 0);
808 }
809
810 if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
811 {
812 *fLog << dbg << "Bin " << i << ": sigmaexc=" << fit.GetSignificanceExc() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
813 fit.PaintResult();
814 }
815 /*
816 if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
817 {
818 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
819 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
820 }*/
821 }
822
823}
824
825
826Bool_t MHAlpha::Finalize()
827{
828 //TH1D *h = fHist.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
829 //h->SetDirectory(0);
830 //Bool_t rc = fFit.Fit(*h);
831 //delete h;
832
833 if (!fFit.FitAlpha(fHist, fOffData))
834 {
835 *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
836 return kTRUE;
837 }
838
839 // Store the final result in fFit
840 *fLog << all;
841 fFit.Print("result");
842
843 fResult->SetVal(fFit.GetMinimizationValue());
844 fSigma->SetVal(fFit.GetGausSigma());
845
846 if (!fSkipHistEnergy)
847 {
848 *fLog << inf << "Processing energy bins..." << endl;
849 FitEnergyBins();
850 }
851 if (!fSkipHistTheta)
852 {
853 *fLog << inf << "Processing theta bins..." << endl;
854 FitThetaBins();
855 }
856 if (!fSkipHistTime)
857 {
858 *fLog << inf << "Processing time bins..." << endl;
859 UpdateAlphaTime(kTRUE);
860 MH::RemoveFirstBin(fHTime);
861 }
862
863 return kTRUE;
864}
865
866// --------------------------------------------------------------------------
867//
868// You can use this function if you want to use a MHMatrix instead of
869// MMcEvt. This function adds all necessary columns to the
870// given matrix. Afterward you should fill the matrix with the corresponding
871// data (eg from a file by using MHMatrix::Fill). If you now loop
872// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
873// will take the values from the matrix instead of the containers.
874//
875// It takes fSkipHist* into account!
876//
877void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
878{
879 if (fMatrix)
880 return;
881
882 fMatrix = mat;
883
884 fMap[0] = fMatrix->AddColumn(GetParameterRule());
885 fMap[1] = -1;
886 fMap[2] = -1;
887 fMap[3] = -1;
888 fMap[4] = -1;
889
890 if (!fSkipHistEnergy)
891 if (type==0)
892 {
893 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
894 fMap[2] = -1;
895 }
896 else
897 {
898 fMap[1] = -1;
899 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
900 }
901
902 if (!fSkipHistTheta)
903 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
904
905 // if (!fSkipHistTime)
906 // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
907}
908
909void MHAlpha::StopMapping()
910{
911 fMatrix = NULL;
912}
913
914Int_t MHAlpha::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
915{
916 Bool_t rc = kFALSE;
917 if (IsEnvDefined(env, prefix, "NumTimeBins", print))
918 {
919 SetNumTimeBins(GetEnvValue(env, prefix, "NumTimeBins", fNumTimeBins));
920 rc = kTRUE;
921 }
922 if (IsEnvDefined(env, prefix, "ForceUsingSize", print))
923 {
924 fForceUsingSize = GetEnvValue(env, prefix, "ForceUsingSize", fForceUsingSize);
925 rc = kTRUE;
926 }
927 return rc;
928}
929
930Int_t MHAlpha::DistancetoPrimitive(Int_t px, Int_t py)
931{
932 // If pad has no subpad return (we are in one of the subpads)
933 if (gPad->GetPad(1)==NULL)
934 return 9999;
935
936 // If pad has a subpad we are in the main pad. Check its value.
937 return gPad->GetPad(1)->DistancetoPrimitive(px,py)==0 ? 0 : 9999;
938}
Note: See TracBrowser for help on using the repository browser.