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

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