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

Last change on this file since 5365 was 5300, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 14.8 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-2004
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#include "MHAlpha.h"
40
41#include <TF1.h>
42#include <TGraph.h>
43#include <TStyle.h>
44#include <TLatex.h>
45#include <TCanvas.h>
46#include <TPaveStats.h>
47#include <TStopwatch.h>
48
49#include "MGeomCam.h"
50#include "MSrcPosCam.h"
51#include "MHillasSrc.h"
52#include "MEnergyEst.h"
53#include "MTime.h"
54#include "MObservatory.h"
55#include "MPointingPos.h"
56#include "MAstroSky2Local.h"
57#include "MStatusDisplay.h"
58#include "MParameters.h"
59#include "MHMatrix.h"
60
61#include "MBinning.h"
62#include "MParList.h"
63
64#include "MLog.h"
65#include "MLogManip.h"
66
67ClassImp(MHAlpha);
68ClassImp(MAlphaFitter);
69
70using namespace std;
71
72// --------------------------------------------------------------------------
73//
74// Default Constructor
75//
76MHAlpha::MHAlpha(const char *name, const char *title)
77 : fResult(0), /*fExcess(0),*/ fEnergy(0), fPointPos(0), fTimeEffOn(0),
78 fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0)
79{
80 //
81 // set the name and title of this object
82 //
83 fName = name ? name : "MHAlpha";
84 fTitle = title ? title : "Alpha plot";
85
86 fHAlpha.SetName("Alpha");
87 fHAlpha.SetTitle("Alpha");
88 fHAlpha.SetXTitle("\\Theta [deg]");
89 fHAlpha.SetYTitle("E_{est} [GeV]");
90 fHAlpha.SetZTitle("|\\alpha| [\\circ]");
91 fHAlpha.SetDirectory(NULL);
92 fHAlpha.UseCurrentStyle();
93
94 // Main histogram
95 fHAlphaTime.SetName("Alpha");
96 fHAlphaTime.SetXTitle("|\\alpha| [\\circ]");
97 fHAlphaTime.SetYTitle("Counts");
98 fHAlphaTime.UseCurrentStyle();
99 fHAlphaTime.SetDirectory(NULL);
100
101
102 fHEnergy.SetName("Energy");
103 fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
104 fHEnergy.SetXTitle("E_{est} [GeV]");
105 fHEnergy.SetYTitle("N_{exc}");
106 fHEnergy.SetDirectory(NULL);
107 fHEnergy.UseCurrentStyle();
108
109 fHTheta.SetName("Theta");
110 fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
111 fHTheta.SetXTitle("\\Theta [\\circ]");
112 fHTheta.SetYTitle("N_{exc}");
113 fHTheta.SetDirectory(NULL);
114 fHTheta.UseCurrentStyle();
115
116 // effective on time versus time
117 fHTime.SetName("Time");
118 fHTime.SetTitle(" N_{exc} vs. Time ");
119 fHTime.SetXTitle("Time");
120 fHTime.SetYTitle("N_{exc} [s]");
121 fHTime.UseCurrentStyle();
122 fHTime.SetDirectory(NULL);
123 fHTime.GetYaxis()->SetTitleOffset(1.2);
124 fHTime.GetXaxis()->SetLabelSize(0.033);
125 fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
126 fHTime.GetXaxis()->SetTimeDisplay(1);
127 fHTime.Sumw2();
128
129 MBinning binsa, binse, binst;
130 binsa.SetEdges(18, 0, 90);
131 binse.SetEdgesLog(25, 10, 100000);
132 binst.SetEdgesCos(50, 0, 60);
133 binse.Apply(fHEnergy);
134 binst.Apply(fHTheta);
135 binsa.Apply(fHAlphaTime);
136
137 MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
138}
139
140void MHAlpha::FitEnergySpec(Bool_t paint)
141{
142 /*
143 TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
144 f.SetParameter(0, -2.1);
145 f.SetParameter(1, 1400); // 50
146 f.SetParameter(2, fHEnergy.Integral()); // 50
147 f.SetLineWidth(2);
148
149 fHEnergy.Fit(&f, "0QI", "", 100, 2400); // Integral? 90, 2800
150
151 const Double_t chi2 = f.GetChisquare();
152 const Int_t NDF = f.GetNDF();
153
154 // was fit successful ?
155 const Bool_t ok = NDF>0 && chi2<2.5*NDF;
156 f.SetLineColor(ok ? kGreen : kRed);
157
158 if (paint)
159 {
160 f.Paint("same");
161
162 if (ok)
163 {
164 TString s;
165 s += Form("\\alpha=%.1f ", f.GetParameter(0));
166 s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
167 s += Form("N=%.1e", f.GetParameter(2));
168
169 TLatex text(0.25, 0.94, s.Data());
170 text.SetBit(TLatex::kTextNDC);
171 text.SetTextSize(0.04);
172 text.Paint();
173 }
174 }*/
175}
176
177void MHAlpha::FitEnergyBins(Bool_t paint)
178{
179 // Do not store the 'final' result in fFit
180 MAlphaFitter fit(fFit);
181
182 const Int_t n = fHAlpha.GetNbinsY();
183
184 for (int i=1; i<=n; i++)
185 {
186 TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":"");
187 if (fit.Fit(*h))
188 {
189 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
190 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
191 }
192 delete h;
193 }
194}
195
196void MHAlpha::FitThetaBins(Bool_t paint)
197{
198 // Do not store the 'final' result in fFit
199 MAlphaFitter fit(fFit);
200
201 const Int_t n = fHAlpha.GetNbinsX();
202
203 for (int i=1; i<=n; i++)
204 {
205 TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":"");
206 if (fit.Fit(*h))
207 {
208 fHTheta.SetBinContent(i, fit.GetEventsExcess());
209 fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
210 }
211 delete h;
212 }
213}
214
215Bool_t MHAlpha::SetupFill(const MParList *pl)
216{
217 fHAlpha.Reset();
218 fHEnergy.Reset();
219 fHTheta.Reset();
220
221 fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
222 if (!fEnergy)
223 *fLog << warn << "MEnergyEst not found... ignored." << endl;
224
225 fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
226 if (!fPointPos)
227 *fLog << warn << "MPointingPos not found... ignored." << endl;
228
229 fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
230 if (!fTimeEffOn)
231 *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
232
233 fTime = (MTime*)pl->FindObject("MTime");
234 if (!fTime)
235 *fLog << warn << "MTime not found... ignored." << endl;
236
237 fResult = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "Significance");
238 if (!fResult)
239 return kFALSE;
240
241 //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
242 //if (!fExcess)
243 // return kFALSE;
244
245 fLastTime = MTime();
246
247 MBinning binst, binse, binsa;
248 binst.SetEdges(fHAlpha, 'x');
249 binse.SetEdges(fHAlpha, 'y');
250 binsa.SetEdges(fHAlpha, 'z');
251
252 MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
253 if (fPointPos && bins)
254 binst.SetEdges(*bins);
255 if (!fPointPos)
256 binst.SetEdges(1, 0, 90);
257
258 bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
259 if (fEnergy && bins)
260 binse.SetEdges(*bins);
261 if (!fEnergy)
262 binse.SetEdges(1, 10, 100000);
263
264 bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
265 if (bins)
266 binsa.SetEdges(*bins);
267
268 binse.Apply(fHEnergy);
269 binst.Apply(fHTheta);
270 binsa.Apply(fHAlphaTime);
271 MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
272
273 MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
274 if (!fit)
275 *fLog << warn << "MAlphaFitter not found... ignored." << endl;
276 else
277 fit->Copy(fFit);
278
279 fFit.Print();
280
281 return kTRUE;
282}
283
284void MHAlpha::InitAlphaTime(const MTime &t)
285{
286 //
287 // If this is the first call we have to initialize the time-histogram
288 //
289 MBinning bins;
290 bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
291 bins.Apply(fHTime);
292
293 fLastTime=t;
294}
295
296void MHAlpha::UpdateAlphaTime(Bool_t final)
297{
298 if (!fTimeEffOn)
299 return;
300
301 const Int_t steps = 10;
302
303 static int rebin = steps;
304
305 if (!final)
306 {
307 if (fLastTime==MTime() && *fTimeEffOn!=MTime())
308 {
309 InitAlphaTime(*fTimeEffOn);
310 return;
311 }
312
313 if (fLastTime!=*fTimeEffOn)
314 {
315 fLastTime=*fTimeEffOn;
316 rebin--;
317 }
318
319 if (rebin>0)
320 return;
321 }
322
323 MAlphaFitter fit(fFit);
324 if (!fit.Fit(fHAlphaTime))
325 return;
326
327 // Reset Histogram
328 fHAlphaTime.Reset();
329
330 // Get number of bins
331 const Int_t n = fHTime.GetNbinsX();
332
333 //
334 // Prepare Histogram
335 //
336
337 MTime *time = final ? fTime : fTimeEffOn;
338 if (final)
339 time->Plus1ns();
340
341 // Enhance binning
342 MBinning bins;
343 bins.SetEdges(fHTime, 'x');
344 bins.AddEdge(time->GetAxisTime());
345 bins.Apply(fHTime);
346
347 //
348 // Fill histogram
349 //
350 fHTime.SetBinContent(n+1, fit.GetEventsExcess());
351 fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);
352
353 *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl;
354
355 rebin = steps;
356}
357
358// --------------------------------------------------------------------------
359//
360// Fill the histogram. For details see the code or the class description
361//
362Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
363{
364 Double_t alpha, energy, theta;
365
366 if (fMatrix)
367 {
368 alpha = (*fMatrix)[fMap[0]];
369 energy = 1000;
370 theta = 0;
371 }
372 else
373 {
374 const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
375 if (!par)
376 {
377 *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
378 return kFALSE;
379 }
380
381 alpha = hil->GetAlpha();
382 energy = fEnergy ? fEnergy->GetEnergy() : 1000;
383 theta = fPointPos ? fPointPos->GetZd() : 0;
384 }
385
386 UpdateAlphaTime();
387
388 fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
389 fHAlphaTime.Fill(TMath::Abs(alpha), w);
390
391 return kTRUE;
392}
393
394// --------------------------------------------------------------------------
395//
396// Paint the integral and the error on top of the histogram
397//
398void MHAlpha::PaintText(Double_t val, Double_t error) const
399{
400 TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
401 text.SetBit(TLatex::kTextNDC);
402 text.SetTextSize(0.04);
403 text.Paint();
404}
405
406// --------------------------------------------------------------------------
407//
408// Update the projections
409//
410void MHAlpha::Paint(Option_t *opt)
411{
412 // Note: Any object cannot be added to a pad twice!
413 // The object is searched by gROOT->FindObject only in
414 // gPad itself!
415 TVirtualPad *padsave = gPad;
416
417 TH1D *h0=0;
418
419 TString o(opt);
420 if (o==(TString)"proj")
421 {
422 TPaveStats *st=0;
423 for (int x=0; x<4; x++)
424 {
425 TVirtualPad *p=padsave->GetPad(x+1);
426 if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
427 continue;
428
429 if (st->GetOptStat()==11)
430 continue;
431
432 const Double_t y1 = st->GetY1NDC();
433 const Double_t y2 = st->GetY2NDC();
434 const Double_t x1 = st->GetX1NDC();
435 const Double_t x2 = st->GetX2NDC();
436
437 st->SetY1NDC((y2-y1)/3+y1);
438 st->SetX1NDC((x2-x1)/3+x1);
439 st->SetOptStat(11);
440 }
441
442 padsave->cd(1);
443 fHAlpha.ProjectionZ(fNameProjAlpha);
444 FitEnergyBins();
445 FitThetaBins();
446 }
447
448 if (o==(TString)"alpha")
449 if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha)))
450 {
451 // Do not store the 'final' result in fFit
452 MAlphaFitter fit(fFit);
453 fit.Fit(*h0, kTRUE);
454 fit.PaintResult();
455 }
456
457 if (o==(TString)"time")
458 PaintText(fHTime.Integral(), 0);
459
460 if (o==(TString)"theta")
461 PaintText(fHTheta.Integral(), 0);
462
463 if (o==(TString)"energy")
464 {
465 if (fHEnergy.GetMaximum()>1)
466 {
467 gPad->SetLogx();
468 gPad->SetLogy();
469 }
470 FitEnergySpec(kTRUE);
471
472 }
473
474 gPad = padsave;
475}
476
477// --------------------------------------------------------------------------
478//
479// Draw the histogram
480//
481void MHAlpha::Draw(Option_t *opt)
482{
483 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
484
485 // Do the projection before painting the histograms into
486 // the individual pads
487 AppendPad("proj");
488
489 pad->SetBorderMode(0);
490 pad->Divide(2,2);
491
492 TH1D *h=0;
493
494 pad->cd(1);
495 gPad->SetBorderMode(0);
496 h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E");
497 h->SetBit(TH1::kNoTitle);
498 h->SetXTitle("\\alpha [\\circ]");
499 h->SetYTitle("Counts");
500 h->SetDirectory(NULL);
501 h->SetMarkerStyle(kFullDotMedium);
502 h->SetBit(kCanDelete);
503 h->Draw();
504 // After the Alpha-projection has been drawn. Fit the histogram
505 // and paint the result into this pad
506 AppendPad("alpha");
507
508 if (fHEnergy.GetNbinsX()>1)
509 {
510 pad->cd(2);
511 gPad->SetBorderMode(0);
512 fHEnergy.Draw();
513 AppendPad("energy");
514 }
515 else
516 delete pad->GetPad(2);
517
518 if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
519 {
520 pad->cd(3);
521 gPad->SetBorderMode(0);
522 fHTime.Draw();
523 AppendPad("time");
524 }
525 else
526 delete pad->GetPad(3);
527
528 if (fHTheta.GetNbinsX()>1)
529 {
530 pad->cd(4);
531 gPad->SetBorderMode(0);
532 fHTheta.Draw();
533 AppendPad("theta");
534 }
535 else
536 delete pad->GetPad(4);
537
538}
539
540Bool_t MHAlpha::Finalize()
541{
542 // Store the final result in fFit
543 TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
544 h->SetDirectory(0);
545 fFit.Print();
546 Bool_t rc = fFit.Fit(*h);
547 delete h;
548 if (!rc)
549 {
550 *fLog << warn << "Histogram empty." << endl;
551 return kTRUE;
552 }
553
554 if (fResult)
555 fResult->SetVal(fFit.GetSignificance());
556
557 FitEnergyBins();
558 FitThetaBins();
559 UpdateAlphaTime(kTRUE);
560 MH::RemoveFirstBin(fHTime);
561
562 return kTRUE;
563}
564
565// --------------------------------------------------------------------------
566//
567// You can use this function if you want to use a MHMatrix instead of
568// MMcEvt. This function adds all necessary columns to the
569// given matrix. Afterward you should fill the matrix with the corresponding
570// data (eg from a file by using MHMatrix::Fill). If you now loop
571// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
572// will take the values from the matrix instead of the containers.
573//
574void MHAlpha::InitMapping(MHMatrix *mat)
575{
576 if (fMatrix)
577 return;
578
579 fMatrix = mat;
580
581 fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
582 //fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
583}
584
585void MHAlpha::StopMapping()
586{
587 fMatrix = NULL;
588}
Note: See TracBrowser for help on using the repository browser.