source: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc@ 7492

Last change on this file since 7492 was 7488, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 33.0 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, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJSpectrum
28//
29// Program to calculate spectrum
30//
31/////////////////////////////////////////////////////////////////////////////
32#include "MJSpectrum.h"
33
34// Root
35#include <TF1.h>
36#include <TH1.h>
37#include <TH2.h>
38#include <TFile.h>
39#include <TChain.h>
40#include <TLatex.h>
41#include <TCanvas.h>
42#include <TObjArray.h>
43
44// Environment
45#include "MLog.h"
46#include "MLogManip.h"
47
48#include "MStatusArray.h"
49#include "MStatusDisplay.h"
50
51// Container
52#include "MH3.h"
53#include "MBinning.h"
54#include "MDataSet.h"
55#include "MMcCorsikaRunHeader.h"
56
57// Spectrum
58#include "../mhflux/MAlphaFitter.h"
59#include "../mhflux/MHAlpha.h"
60#include "../mhflux/MHCollectionArea.h"
61//#include "../mhflux/MHThreshold.h"
62#include "../mhflux/MHEnergyEst.h"
63#include "../mhflux/MMcSpectrumWeight.h"
64
65// Eventloop
66#include "MEvtLoop.h"
67#include "MTaskList.h"
68#include "MParList.h"
69
70// Tasks/Filter
71#include "MReadMarsFile.h"
72#include "MReadMarsFile.h"
73#include "MFEventSelector2.h"
74#include "MFDataMember.h"
75#include "MEnergyEstimate.h"
76#include "MTaskEnv.h"
77#include "MFillH.h"
78#include "MHillasCalc.h"
79//#include "MSrcPosCalc.h"
80#include "MContinue.h"
81
82ClassImp(MJSpectrum);
83
84using namespace std;
85
86MJSpectrum::MJSpectrum(const char *name, const char *title)
87 : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
88 fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE),
89 fNoThetaWeights(kFALSE)
90{
91 fName = name ? name : "MJSpectrum";
92 fTitle = title ? title : "Standard program to calculate spectrum";
93}
94
95MJSpectrum::~MJSpectrum()
96{
97 if (fCut0)
98 delete fCut0;
99 if (fCut1)
100 delete fCut1;
101 if (fCut2)
102 delete fCut2;
103 if (fCut3)
104 delete fCut3;
105 if (fEstimateEnergy)
106 delete fEstimateEnergy;
107}
108
109// --------------------------------------------------------------------------
110//
111// Setup a task estimating the energy. The given task is cloned.
112//
113void MJSpectrum::SetEnergyEstimator(const MTask *task)
114{
115 if (fEstimateEnergy)
116 delete fEstimateEnergy;
117 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
118}
119
120Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
121{
122 if (task)
123 {
124 delete task;
125 task = 0;
126 }
127
128 task = (MTask*)gFile->Get(name);
129 if (!task)
130 {
131 *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
132 return kFALSE;
133 }
134 if (!task->InheritsFrom(MTask::Class()))
135 {
136 *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
137 delete task;
138 return kFALSE;
139 }
140
141 task->SetName(name);
142
143 if (dynamic_cast<MContinue*>(task))
144 dynamic_cast<MContinue*>(task)->SetAllowEmpty();
145
146 return kTRUE;
147}
148
149void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
150{
151 fLog->Separator("Alpha Fitter");
152 *fLog << all;
153 fit.Print();
154
155 fLog->Separator("Used Cuts");
156 fCut0->Print();
157 fCut1->Print();
158 fCut2->Print();
159 fCut3->Print();
160
161 //gLog.Separator("Energy Estimator");
162 //fEstimateEnergy->Print();
163}
164
165// --------------------------------------------------------------------------
166//
167// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
168// the dataset.
169// The simulated energy range and spectral slope is initialized from
170// there.
171// In the following eventloops the forced check in MMcSpectrumWeight
172// ensures, that the spectral slope and energy range doesn't change.
173//
174Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
175{
176 fLog->Separator("Initialize energy weighting");
177
178 if (!CheckEnv(w))
179 {
180 *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
181 return kFALSE;
182 }
183
184 TChain chain("RunHeaders");
185 set.AddFilesOn(chain);
186
187 MMcCorsikaRunHeader *h=0;
188 chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
189 chain.GetEntry(1);
190
191 if (!h)
192 {
193 *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
194 return kFALSE;
195 }
196
197 if (!w.Set(*h))
198 {
199 *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
200 return kFALSE;
201 }
202
203 w.Print();
204 return kTRUE;
205}
206
207Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
208{
209 *fLog << inf << "Reading from file: " << fPathIn << endl;
210
211 TFile file(fPathIn, "READ");
212 if (!file.IsOpen())
213 {
214 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
215 return -1;
216 }
217
218 MStatusArray arr;
219 if (arr.Read()<=0)
220 {
221 *fLog << "MStatusDisplay not found in file... abort." << endl;
222 return -1;
223 }
224
225 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
226 TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
227 if (!vstime || !size)
228 return -1;
229
230 vstime->Copy(h1);
231 size->Copy(h2);
232 h1.SetDirectory(0);
233 h2.SetDirectory(0);
234
235 if (fDisplay)
236 arr.DisplayIn(*fDisplay, "Hist");
237
238 if (!ReadTask(fCut0, "Cut0"))
239 return -1;
240 if (!ReadTask(fCut1, "Cut1"))
241 return -1;
242 if (!ReadTask(fCut2, "Cut2"))
243 return -1;
244 if (!ReadTask(fCut3, "Cut3"))
245 return -1;
246
247 TObjArray arrread;
248
249 TIter Next(plist);
250 TObject *o=0;
251 while ((o=Next()))
252 if (o->InheritsFrom(MBinning::Class()))
253 arrread.Add(o);
254
255 arrread.Add(plist.FindObject("MAlphaFitter"));
256
257 if (!ReadContainer(arrread))
258 return -1;
259
260 return vstime->Integral();
261}
262
263Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
264{
265 // Some debug output
266 fLog->Separator("Compiling original MC distribution");
267
268 weight.SetNameMcEvt("MMcEvtBasic");
269 const TString w(weight.GetFormulaWeights());
270 weight.SetNameMcEvt();
271
272 *fLog << inf << "Using weights: " << w << endl;
273 *fLog << "Please stand by, this may take a while..." << flush;
274
275 if (fDisplay)
276 fDisplay->SetStatusLine1("Compiling MC distribution...");
277
278 // Create chain
279 TChain chain("OriginalMC");
280 set.AddFilesOn(chain);
281
282 // Prepare histogram
283 h.Reset();
284
285 // Fill histogram from chain
286 h.SetDirectory(gROOT);
287 if (h.InheritsFrom(TH2::Class()))
288 {
289 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
290 h.SetXTitle("\\Theta [\\circ]");
291 h.SetYTitle("E [GeV]");
292 h.SetZTitle("Counts");
293 chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
294 }
295 else
296 {
297 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
298 h.SetXTitle("\\Theta [\\circ]");
299 h.SetYTitle("Counts");
300 chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
301 }
302 h.SetDirectory(0);
303
304 *fLog << "done." << endl;
305 if (fDisplay)
306 fDisplay->SetStatusLine2("done.");
307
308 if (h.GetEntries()>0)
309 return kTRUE;
310
311 *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
312
313 return h.GetEntries()>0;
314}
315
316Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
317{
318 // Display some stuff
319 if (fDisplay)
320 {
321 TCanvas &c = fDisplay->AddTab("ZdDist");
322 c.Divide(2,2);
323
324 // On-Time vs. Theta
325 c.cd(1);
326 gPad->SetBorderMode(0);
327 temp1.DrawCopy();
328
329 // Number of MC events (produced) vs Theta
330 c.cd(2);
331 gPad->SetBorderMode(0);
332 temp2.SetName("NVsTheta");
333 temp2.DrawCopy();
334
335 c.cd(4);
336 gPad->SetBorderMode(0);
337
338 c.cd(3);
339 gPad->SetBorderMode(0);
340 }
341
342 // Calculate the Probability
343 temp1.Divide(&temp2);
344 temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());
345
346 // Some cosmetics: Name, Axis, etc.
347 temp1.SetName("ProbVsTheta");
348 temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
349 temp1.SetYTitle("Probability");
350 if (fDisplay)
351 temp1.DrawCopy();
352
353 return kTRUE;
354}
355
356// --------------------------------------------------------------------------
357//
358// Display the final theta distribution.
359//
360void MJSpectrum::DisplayResult(const TH2D &h2) const
361{
362 if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
363 return;
364
365 TH1D &proj = *h2.ProjectionX();
366 proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
367 proj.SetXTitle("\\Theta [\\circ]");
368 proj.SetYTitle("Counts");
369 proj.SetLineColor(kBlue);
370 proj.SetDirectory(0);
371 proj.SetBit(kCanDelete);
372
373 TVirtualPad *pad = gPad;
374
375 pad->cd(4);
376 proj.DrawCopy();
377
378 pad->cd(1);
379 TH1D *theta = (TH1D*)gPad->FindObject("Theta");
380 if (theta)
381 {
382 proj.Scale(theta->GetMaximum()/proj.GetMaximum());
383 theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
384 }
385 proj.Draw("same");
386}
387
388// --------------------------------------------------------------------------
389//
390// Fills the excess histogram (vs E-est) from the events stored in the
391// ganymed result file and therefor estimates the energy.
392//
393// The resulting histogram excess-vs-energy ist copied into h2.
394//
395Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
396{
397 // Try to find the class used to determin the signal!
398 TString cls("MHAlpha");
399 if (fDisplay)
400 {
401 TCanvas *c = fDisplay->GetCanvas("Hist");
402 if (c)
403 {
404 TIter Next(c->GetListOfPrimitives());
405 TObject *obj=0;
406 while ((obj=Next()))
407 if (obj->InheritsFrom(MHAlpha::Class()))
408 break;
409 if (obj)
410 cls = obj->ClassName();
411 }
412 }
413
414 cout << "FOUND: "<< cls << endl;
415
416 // Now fill the histogram
417 *fLog << endl;
418 fLog->Separator("Refill Excess");
419 *fLog << endl;
420
421 MTaskList tlist;
422 plist.AddToList(&tlist);
423
424 MReadTree read("Events");
425 read.DisableAutoScheme();
426 read.AddFile(fPathIn);
427
428 MEnergyEstimate est;
429 MTaskEnv taskenv1("EstimateEnergy");
430 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
431
432 // FIXME: Create HistE and HistEOff to be able to modify it from
433 // the resource file.
434
435 MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
436 MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
437
438 MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
439 MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
440
441 fill1.SetFilter(&f0);
442 fill2.SetFilter(&f1);
443
444 tlist.AddToList(&read);
445 tlist.AddToList(&taskenv1);
446 tlist.AddToList(&f0);
447 tlist.AddToList(&f1);
448 tlist.AddToList(&fill1);
449 tlist.AddToList(&fill2);
450
451 MEvtLoop loop("RefillExcess"); // ***** fName *****
452 loop.SetParList(&plist);
453 loop.SetDisplay(fDisplay);
454 loop.SetLogStream(fLog);
455
456 if (!SetupEnv(loop))
457 return kFALSE;
458
459 if (!loop.Eventloop())
460 {
461 *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
462 return kFALSE;
463 }
464
465 if (!loop.GetDisplay())
466 {
467 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
468 return kFALSE;
469 }
470
471 const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
472 if (!halpha)
473 {
474 *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
475 return kFALSE;
476 }
477
478 halpha->GetHEnergy().Copy(h2);
479 h2.SetDirectory(0);
480
481 return kTRUE;
482}
483
484Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
485{
486 MTaskList tlist1;
487 plist.Replace(&tlist1);
488
489 MReadMarsFile readmc("OriginalMC");
490 //readmc.DisableAutoScheme();
491 set.AddFilesOn(readmc);
492 readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
493 readmc.EnableBranch("MMcEvtBasic.fEnergy");
494
495 mh1.SetLogy();
496 mh1.SetLogz();
497 mh1.SetName("ThetaE");
498
499 MFillH fill0(&mh1);
500 //fill0.SetDrawOption("projx only");
501
502 MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
503 MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
504 if (bins2 && bins3)
505 {
506 bins2->SetName("BinningThetaEY");
507 bins3->SetName("BinningThetaEX");
508 }
509 tlist1.AddToList(&readmc);
510 tlist1.AddToList(&weight);
511
512 temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
513 MH3 mh3mc(temp1);
514
515 MFEventSelector2 sel1(mh3mc);
516 sel1.SetHistIsProbability();
517
518 fill0.SetFilter(&sel1);
519
520 if (!fRawMc)
521 tlist1.AddToList(&sel1);
522 tlist1.AddToList(&fill0);
523
524 MEvtLoop loop1("IntermediateLoop"); // ***** fName *****
525 loop1.SetParList(&plist);
526 loop1.SetLogStream(fLog);
527 loop1.SetDisplay(fDisplay);
528
529 if (!SetupEnv(loop1))
530 return kFALSE;
531
532 if (!loop1.Eventloop(fMaxEvents))
533 {
534 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
535 return kFALSE;
536 }
537
538 if (!loop1.GetDisplay())
539 {
540 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
541 return kFALSE;
542 }
543
544 if (bins2 && bins3)
545 {
546 bins2->SetName("BinningEnergyEst");
547 bins3->SetName("BinningTheta");
548 }
549
550 return kTRUE;
551}
552
553TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
554{
555 const Double_t p0 = f.GetParameter(0);
556 const Double_t p1 = f.GetParameter(1);
557
558 const Double_t e0 = f.GetParError(0);
559 const Double_t e1 = f.GetParError(1);
560
561 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
562 const Double_t exp = TMath::Power(10, np);
563
564 const Float_t fe0 = TMath::Log10(e0);
565 const Float_t fe1 = TMath::Log10(e1);
566
567 // calc number of (gueltige ziffern) digits to be displayed
568 const char *f0 = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0";
569 const char *f1 = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0";
570
571 TString str;
572 switch (type)
573 {
574 case 0:
575 {
576 const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1);
577 const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
578
579 str = Form(fmt0.Data(), p1/exp, e1/exp, np);
580 str += Form(fmt1.Data(), p0, e0);
581 str += "\\frac{ph}{TeVm^{2}s}";
582 }
583 break;
584 case 1:
585 str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
586 break;
587 case 2:
588 str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF()));
589 break;
590 }
591 return str;
592}
593
594TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
595{
596 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
597 f.SetParameter(0, -2.5);
598 f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000)));
599 f.SetLineColor(kBlue);
600 spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped
601 f.DrawCopy("same");
602
603 TString str = FormString(f);
604
605 TLatex tex;
606 tex.SetTextSize(0.045);
607 tex.SetBit(TLatex::kTextNDC);
608 tex.SetTextAlign(31);
609 tex.DrawLatex(0.89, 0.935, str);
610
611 str = FormString(f, 1);
612 tex.DrawLatex(0.89, 0.83, str);
613
614 str = FormString(f, 2);
615 tex.DrawLatex(0.89, 0.735, str);
616
617 TArrayD res(2);
618 res[0] = f.GetParameter(0);
619 res[1] = f.GetParameter(1);
620 return res;
621}
622
623// --------------------------------------------------------------------------
624//
625// Calculate the final spectrum from:
626// - collection area
627// - excess
628// - correction coefficients
629// - ontime
630// and display it
631//
632TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
633{
634 TH1D collarea(area.GetHEnergy());
635 TH1D spectrum(excess);
636 TH1D weights;
637 hest.GetWeights(weights);
638
639 cout << "Effective On time: " << ontime << "s" << endl;
640
641 spectrum.SetDirectory(NULL);
642 spectrum.SetBit(kCanDelete);
643 spectrum.Scale(1./ontime);
644 spectrum.Divide(&collarea);
645 spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
646 spectrum.SetYTitle("N/sm^{2}");
647
648 TCanvas &c1 = fDisplay->AddTab("Spectrum");
649 c1.Divide(2,2);
650 c1.cd(1);
651 gPad->SetBorderMode(0);
652 gPad->SetLogx();
653 gPad->SetLogy();
654 gPad->SetGridx();
655 gPad->SetGridy();
656 collarea.DrawCopy();
657
658 c1.cd(2);
659 gPad->SetBorderMode(0);
660 gPad->SetLogx();
661 gPad->SetLogy();
662 gPad->SetGridx();
663 gPad->SetGridy();
664 TH1D *spec=(TH1D*)spectrum.DrawCopy();
665 FitSpectrum(*spec);
666
667 c1.cd(3);
668 gPad->SetBorderMode(0);
669 gPad->SetLogx();
670 gPad->SetLogy();
671 gPad->SetGridx();
672 gPad->SetGridy();
673 weights.DrawCopy();
674
675 //spectrum.Multiply(&weights);
676 spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
677 spectrum.SetBit(TH1::kNoStats);
678
679 for (int i=0; i<excess.GetNbinsX(); i++)
680 {
681 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
682 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
683
684 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
685 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
686 }
687
688 c1.cd(4);
689 gPad->SetBorderMode(0);
690 gPad->SetLogx();
691 gPad->SetLogy();
692 gPad->SetGridx();
693 gPad->SetGridy();
694 spectrum.SetMinimum(1e-12);
695 spectrum.SetXTitle("E [GeV]");
696 spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
697 spec = (TH1D*)spectrum.DrawCopy();
698
699 // Calculate Spectrum * E^2
700 for (int i=0; i<spectrum.GetNbinsX(); i++)
701 {
702 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
703 spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) *spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
704 }
705
706 // Display dN/dE*E^2 for conveinience
707 fDisplay->AddTab("Check");
708 gPad->SetLogx();
709 gPad->SetLogy();
710 gPad->SetGrid();
711
712 spectrum.SetName("FluxStd");
713 spectrum.SetMarkerStyle(kFullDotMedium);
714 spectrum.SetTitle("Differential flux times E^{2}");
715 spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
716 spectrum.SetDirectory(0);
717 spectrum.DrawCopy();
718
719 // Fit Spectrum
720 c1.cd(4);
721 return FitSpectrum(*spec);
722
723/*
724 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
725 f.SetParameter(0, -2.87);
726 f.SetParameter(1, 1.9e-6);
727 f.SetLineColor(kGreen);
728 spectrum.Fit(&f, "NIM", "", 100, 5000);
729 f.DrawCopy("same");
730
731 const Double_t p0 = f.GetParameter(0);
732 const Double_t p1 = f.GetParameter(1);
733
734 const Double_t e0 = f.GetParError(0);
735 const Double_t e1 = f.GetParError(1);
736
737 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
738 const Double_t exp = TMath::Power(10, np);
739
740 TString str;
741 str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
742 str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
743 str += "\\frac{ph}{TeVm^{2}s}";
744
745 TLatex tex;
746 tex.SetTextSize(0.045);
747 tex.SetBit(TLatex::kTextNDC);
748 tex.DrawLatex(0.45, 0.935, str);
749
750 str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
751 tex.DrawLatex(0.70, 0.83, str);
752
753 TArrayD res(2);
754 res[0] = f.GetParameter(0);
755 res[1] = f.GetParameter(1);
756
757 return res;
758 */
759/*
760 // Plot other spectra from Whipple
761 f.SetParameter(0, -2.45);
762 f.SetParameter(1, 3.3e-7);
763 f.SetRange(300, 8000);
764 f.SetLineColor(kBlack);
765 f.SetLineStyle(kDashed);
766 f.DrawCopy("same");
767
768 // Plot other spectra from Cangaroo
769 f.SetParameter(0, -2.53);
770 f.SetParameter(1, 2.0e-7);
771 f.SetRange(7000, 50000);
772 f.SetLineColor(kBlack);
773 f.SetLineStyle(kDashed);
774 f.DrawCopy("same");
775
776 // Plot other spectra from Robert
777 f.SetParameter(0, -2.59);
778 f.SetParameter(1, 2.58e-7);
779 f.SetRange(150, 1500);
780 f.SetLineColor(kBlack);
781 f.SetLineStyle(kDashed);
782 f.DrawCopy("same");
783
784 // Plot other spectra from HEGRA
785 f.SetParameter(0, -2.61);
786 f.SetParameter(1, 2.7e-7);
787 f.SetRange(1000, 20000);
788 f.SetLineColor(kBlack);
789 f.SetLineStyle(kDashed);
790 f.DrawCopy("same");
791 */
792}
793
794// --------------------------------------------------------------------------
795//
796// Scale some image parameter plots using the scale factor and plot them
797// together with the corresponding MC histograms.
798// Called from DisplaySize
799//
800Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
801{
802 TString same(name);
803 same += "Same";
804
805 TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
806 TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
807 if (!h1 || !h2)
808 return kFALSE;
809
810 TObject *obj = plist.FindObject(plot);
811 if (!obj)
812 {
813 *fLog << warn << plot << " not in parameter list... skipping." << endl;
814 return kFALSE;
815 }
816
817 TH1 *h3 = (TH1*)obj->FindObject(name);
818 if (!h3)
819 {
820 *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
821 return kFALSE;
822 }
823
824
825 const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
826 const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
827
828 gPad->SetBorderMode(0);
829 h2->SetLineColor(kBlack);
830 h3->SetLineColor(kBlue);
831 h2->Add(h1, -ascale);
832
833 //h2->Scale(1./ontime); //h2->Integral());
834 h3->Scale(scale); //h3->Integral());
835
836 h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
837
838 h2 = h2->DrawCopy();
839 h3 = h3->DrawCopy("same");
840
841 // Don't do this on the original object!
842 h2->SetStats(kFALSE);
843 h3->SetStats(kFALSE);
844
845 return kTRUE;
846}
847
848// --------------------------------------------------------------------------
849//
850// Take a lot of histograms and plot them together in one plot.
851// Calls PlotSame
852//
853Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
854{
855 *fLog << inf << "Reading from file: " << fPathIn << endl;
856
857 TFile file(fPathIn, "READ");
858 if (!file.IsOpen())
859 {
860 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
861 return kFALSE;
862 }
863
864 file.cd();
865 MStatusArray arr;
866 if (arr.Read()<=0)
867 {
868 *fLog << "MStatusDisplay not found in file... abort." << endl;
869 return kFALSE;
870 }
871
872 TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
873 if (!excess)
874 return kFALSE;
875
876 // ------------------- Plot excess versus size -------------------
877
878 TCanvas &c = fDisplay->AddTab("Excess");
879 c.Divide(3,2);
880 c.cd(1);
881 gPad->SetBorderMode(0);
882 gPad->SetLogx();
883 gPad->SetLogy();
884 gPad->SetGridx();
885 gPad->SetGridy();
886
887 excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
888 excess = excess->DrawCopy("E2");
889 // Don't do this on the original object!
890 excess->SetStats(kFALSE);
891 excess->SetMarkerStyle(kFullDotMedium);
892 excess->SetFillColor(kBlack);
893 excess->SetFillStyle(0);
894 excess->SetName("Excess ");
895 excess->SetDirectory(0);
896
897 TObject *o=0;
898 if ((o=plist.FindObject("ExcessMC")))
899 {
900 TH1 *histsel = (TH1F*)o->FindObject("");
901 if (histsel)
902 {
903 if (scale<0)
904 scale = excess->Integral()/histsel->Integral();
905
906 histsel->Scale(scale);
907 histsel->SetLineColor(kBlue);
908 histsel->SetBit(kCanDelete);
909 histsel = histsel->DrawCopy("E1 same");
910 // Don't do this on the original object!
911 histsel->SetStats(kFALSE);
912
913 fLog->Separator("Kolmogorov Test");
914 histsel->KolmogorovTest(excess, "DX");
915 fLog->Separator("Chi^2 Test");
916 const Double_t p = histsel->Chi2Test(excess, "P");
917
918 TLatex tex;
919 tex.SetBit(TLatex::kTextNDC);
920 tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
921 }
922 }
923
924 // -------------- Comparison of Image Parameters --------------
925 c.cd(2);
926 PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
927
928 c.cd(3);
929 PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
930
931 c.cd(4);
932 PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
933
934 c.cd(5);
935 PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
936
937 c.cd(6);
938 PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
939
940 return kTRUE;
941}
942
943Bool_t MJSpectrum::Process(const MDataSet &set)
944{
945 if (!set.IsValid())
946 {
947 *fLog << err << "ERROR - DataSet invalid!" << endl;
948 return kFALSE;
949 }
950
951 CheckEnv();
952
953 // --------------------------------------------------------------------------------
954
955 *fLog << inf;
956 fLog->Separator(GetDescriptor());
957 *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
958 *fLog << endl;
959
960 MBinning bins1("BinningAlpha");
961 MBinning bins2("BinningEnergyEst");
962 MBinning bins3("BinningTheta");
963 MBinning bins4("BinningFalseSource");
964 MBinning bins5("BinningWidth");
965 MBinning bins6("BinningLength");
966 MBinning bins7("BinningDist");
967 MBinning bins8("BinningMaxDist");
968 MBinning bins9("BinningM3Long");
969 MBinning bins0("BinningConc1");
970
971 MAlphaFitter fit;
972
973 MParList plist;
974 plist.AddToList(&bins1);
975 plist.AddToList(&bins3);
976 plist.AddToList(&bins4);
977 plist.AddToList(&bins5);
978 plist.AddToList(&bins6);
979 plist.AddToList(&bins7);
980 plist.AddToList(&bins8);
981 plist.AddToList(&bins9);
982 plist.AddToList(&bins0);
983 plist.AddToList(&fit);
984
985 TH1D temp1, size;
986 const Float_t ontime = ReadInput(plist, temp1, size);
987 if (ontime<0)
988 {
989 *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
990 return kFALSE;
991 }
992
993 plist.AddToList(&bins2);
994
995 MMcSpectrumWeight weight;
996 if (!InitWeighting(set, weight))
997 return kFALSE;
998
999 PrintSetup(fit);
1000 bins3.SetEdges(temp1, 'x');
1001
1002 TH1D temp2(temp1);
1003 if (!ReadOrigMCDistribution(set, temp2, weight))
1004 return kFALSE;
1005
1006 if (!GetThetaDistribution(temp1, temp2))
1007 return kFALSE;
1008
1009 if (!fNoThetaWeights)
1010 weight.SetWeightsZd(&temp1);
1011
1012 TH1D excess;
1013 if (!Refill(plist, excess))
1014 return kFALSE;
1015
1016 TH2D hist;
1017 MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
1018 if (fSimpleMode)
1019 {
1020 hist.UseCurrentStyle();
1021 MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
1022 if (!ReadOrigMCDistribution(set, hist, weight))
1023 return kFALSE;
1024
1025 if (!fRawMc)
1026 {
1027 for (int y=0; y<hist.GetNbinsY(); y++)
1028 for (int x=0; x<hist.GetNbinsX(); x++)
1029 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
1030 //hist.SetEntries(hist.Integral());
1031 }
1032 }
1033 else
1034 {
1035 weight.SetNameMcEvt("MMcEvtBasic");
1036 if (!IntermediateLoop(plist, mh1, temp1, set, weight))
1037 return kFALSE;
1038 weight.SetNameMcEvt();
1039 }
1040
1041 DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1042
1043 // ------------------------- Final loop --------------------------
1044
1045 *fLog << endl;
1046 fLog->Separator("Calculate efficiencies");
1047 *fLog << endl;
1048
1049 MTaskList tlist2;
1050 plist.Replace(&tlist2);
1051
1052 MReadMarsFile read("Events");
1053 read.DisableAutoScheme();
1054 set.AddFilesOn(read);
1055
1056 // Selector to get correct (final) theta-distribution
1057 temp1.SetXTitle("MPointingPos.fZd");
1058
1059 MH3 mh3(temp1);
1060
1061 MFEventSelector2 sel2(mh3);
1062 sel2.SetHistIsProbability();
1063
1064 MContinue contsel(&sel2);
1065 contsel.SetInverted();
1066
1067 // Get correct source position
1068 //MSrcPosCalc calc;
1069
1070 // Calculate corresponding Hillas parameters
1071 /*
1072 MHillasCalc hcalc1;
1073 MHillasCalc hcalc2("MHillasCalcAnti");
1074 hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1075 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1076 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1077 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1078 */
1079
1080 // Fill collection area and energy estimator (unfolding)
1081 // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1082 MHCollectionArea area;
1083 area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1084 MHEnergyEst hest;
1085
1086 MFillH fill3(&area, "", "FillCollectionArea");
1087 MFillH fill4(&hest, "", "FillEnergyEst");
1088 MFillH fill5("MHThreshold", "", "FillThreshold");
1089 fill3.SetWeight();
1090 fill4.SetWeight();
1091 fill5.SetWeight();
1092 fill3.SetNameTab("ColArea");
1093 fill4.SetNameTab("E-Est");
1094 fill5.SetNameTab("Threshold");
1095
1096 MH3 hsize("MHillas.fSize");
1097 hsize.SetName("ExcessMC");
1098 hsize.Sumw2();
1099
1100 MBinning bins(size, "BinningExcessMC");
1101 plist.AddToList(&hsize);
1102 plist.AddToList(&bins);
1103
1104 MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1105 MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1106 MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1107 MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1108 MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1109 MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1110 MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1111 MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1112 fill1a.SetNameTab("PreCut");
1113 fill2a.SetNameTab("PostCut");
1114 fill3a.SetNameTab("VsSize");
1115 fill4a.SetNameTab("HilExt");
1116 fill5a.SetNameTab("HilSrc");
1117 fill6a.SetNameTab("ImgPar");
1118 fill7a.SetNameTab("NewPar");
1119 fill8a.SetBit(MFillH::kDoNotDisplay);
1120 fill1a.SetWeight();
1121 fill2a.SetWeight();
1122 fill3a.SetWeight();
1123 fill4a.SetWeight();
1124 fill5a.SetWeight();
1125 fill6a.SetWeight();
1126 fill7a.SetWeight();
1127 fill8a.SetWeight();
1128
1129 MEnergyEstimate est;
1130 MTaskEnv taskenv1("EstimateEnergy");
1131 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1132
1133 tlist2.AddToList(&read);
1134 if (!fRawMc && fNoThetaWeights)
1135 tlist2.AddToList(&contsel);
1136 //tlist2.AddToList(&calc);
1137 //tlist2.AddToList(&hcalc1);
1138 //tlist2.AddToList(&hcalc2);
1139 tlist2.AddToList(&weight);
1140 tlist2.AddToList(&fill1a);
1141 tlist2.AddToList(fCut0);
1142 tlist2.AddToList(fCut1);
1143 tlist2.AddToList(fCut2);
1144 tlist2.AddToList(fCut3);
1145 tlist2.AddToList(&taskenv1);
1146 tlist2.AddToList(&fill3);
1147 tlist2.AddToList(&fill4);
1148 tlist2.AddToList(&fill5);
1149 tlist2.AddToList(&fill2a);
1150 tlist2.AddToList(&fill3a);
1151 tlist2.AddToList(&fill4a);
1152 tlist2.AddToList(&fill5a);
1153 tlist2.AddToList(&fill6a);
1154 tlist2.AddToList(&fill7a);
1155 tlist2.AddToList(&fill8a);
1156 //tlist2.AddToList(&fill9a);
1157
1158 MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
1159 loop2.SetParList(&plist);
1160 loop2.SetDisplay(fDisplay);
1161 loop2.SetLogStream(fLog);
1162
1163 if (!SetupEnv(loop2))
1164 return kFALSE;
1165
1166 if (!loop2.Eventloop(fMaxEvents))
1167 {
1168 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1169 return kFALSE;
1170 }
1171
1172 if (!loop2.GetDisplay())
1173 {
1174 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1175 return kFALSE;
1176 }
1177
1178 gLog.Separator("Energy Estimator");
1179 if (plist.FindObject("EstimateEnergy"))
1180 plist.FindObject("EstimateEnergy")->Print();
1181
1182 gLog.Separator("Spectrum");
1183
1184 // -------------------------- Spectrum ----------------------------
1185
1186 // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1187 TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
1188
1189 // Spectrum fitted (convert res[1] from TeV to GeV)
1190 TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
1191
1192 // Number of events this spectrum would produce per s and m^2
1193 Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1194
1195 // scale with effective collection area to get the event rate (N/s)
1196 // scale with the effective observation time to absolute observed events
1197 n *= area.GetCollectionAreaAbs()*ontime; // N
1198
1199 // Now calculate the scale factor from the number of events
1200 // produced and the number of events which should have been
1201 // observed with our telescope in the time ontime
1202 const Double_t scale = n/area.GetEntries();
1203
1204 // Print normalization constant
1205 cout << "MC normalization factor: " << scale << endl;
1206
1207 // Overlay normalized plots
1208 DisplaySize(plist, scale);
1209
1210 // check if output should be written
1211 if (!fPathOut.IsNull())
1212 fDisplay->SaveAsRoot(fPathOut);
1213
1214 return kTRUE;
1215}
Note: See TracBrowser for help on using the repository browser.