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

Last change on this file since 8006 was 7650, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 33.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, 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 // by setting it here it is distributed to all consecutive tasks
452 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
453
454 MEvtLoop loop("RefillExcess"); // ***** fName *****
455 loop.SetParList(&plist);
456 loop.SetDisplay(fDisplay);
457 loop.SetLogStream(fLog);
458
459 if (!SetupEnv(loop))
460 return kFALSE;
461
462 if (!loop.Eventloop())
463 {
464 *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
465 return kFALSE;
466 }
467
468 if (!loop.GetDisplay())
469 {
470 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
471 return kFALSE;
472 }
473
474 const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
475 if (!halpha)
476 {
477 *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
478 return kFALSE;
479 }
480
481 halpha->GetHEnergy().Copy(h2);
482 h2.SetDirectory(0);
483
484 return kTRUE;
485}
486
487Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
488{
489 MTaskList tlist1;
490 plist.Replace(&tlist1);
491
492 MReadMarsFile readmc("OriginalMC");
493 //readmc.DisableAutoScheme();
494 set.AddFilesOn(readmc);
495 readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
496 readmc.EnableBranch("MMcEvtBasic.fEnergy");
497
498 mh1.SetLogy();
499 mh1.SetLogz();
500 mh1.SetName("ThetaE");
501
502 MFillH fill0(&mh1);
503 //fill0.SetDrawOption("projx only");
504
505 MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
506 MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
507 if (bins2 && bins3)
508 {
509 bins2->SetName("BinningThetaEY");
510 bins3->SetName("BinningThetaEX");
511 }
512 tlist1.AddToList(&readmc);
513 tlist1.AddToList(&weight);
514
515 temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
516 MH3 mh3mc(temp1);
517
518 MFEventSelector2 sel1(mh3mc);
519 sel1.SetHistIsProbability();
520
521 fill0.SetFilter(&sel1);
522
523 if (!fRawMc)
524 tlist1.AddToList(&sel1);
525 tlist1.AddToList(&fill0);
526
527 // by setting it here it is distributed to all consecutive tasks
528 tlist1.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
529
530 MEvtLoop loop1("IntermediateLoop"); // ***** fName *****
531 loop1.SetParList(&plist);
532 loop1.SetLogStream(fLog);
533 loop1.SetDisplay(fDisplay);
534
535 if (!SetupEnv(loop1))
536 return kFALSE;
537
538 if (!loop1.Eventloop(fMaxEvents))
539 {
540 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
541 return kFALSE;
542 }
543
544 if (!loop1.GetDisplay())
545 {
546 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
547 return kFALSE;
548 }
549
550 if (bins2 && bins3)
551 {
552 bins2->SetName("BinningEnergyEst");
553 bins3->SetName("BinningTheta");
554 }
555
556 return kTRUE;
557}
558
559TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
560{
561 const Double_t p0 = f.GetParameter(0);
562 const Double_t p1 = f.GetParameter(1);
563
564 const Double_t e0 = f.GetParError(0);
565 const Double_t e1 = f.GetParError(1);
566
567 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
568 const Double_t exp = TMath::Power(10, np);
569
570 const Float_t fe0 = TMath::Log10(e0);
571 const Float_t fe1 = TMath::Log10(e1);
572
573 // calc number of (gueltige ziffern) digits to be displayed
574 const char *f0 = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0";
575 const char *f1 = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0";
576
577 TString str;
578 switch (type)
579 {
580 case 0:
581 {
582 const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1);
583 const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
584
585 str = Form(fmt0.Data(), p1/exp, e1/exp, np);
586 str += Form(fmt1.Data(), p0, e0);
587 str += "\\frac{ph}{TeVm^{2}s}";
588 }
589 break;
590 case 1:
591 str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
592 break;
593 case 2:
594 str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF()));
595 break;
596 }
597 return str;
598}
599
600TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
601{
602 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
603 f.SetParameter(0, -2.5);
604 f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000)));
605 f.SetLineColor(kBlue);
606 spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped
607 f.DrawCopy("same");
608
609 TString str = FormString(f);
610
611 TLatex tex;
612 tex.SetTextSize(0.045);
613 tex.SetBit(TLatex::kTextNDC);
614 tex.SetTextAlign(31);
615 tex.DrawLatex(0.89, 0.935, str);
616
617 str = FormString(f, 1);
618 tex.DrawLatex(0.89, 0.83, str);
619
620 str = FormString(f, 2);
621 tex.DrawLatex(0.89, 0.735, str);
622
623 TArrayD res(2);
624 res[0] = f.GetParameter(0);
625 res[1] = f.GetParameter(1);
626 return res;
627}
628
629// --------------------------------------------------------------------------
630//
631// Calculate the final spectrum from:
632// - collection area
633// - excess
634// - correction coefficients
635// - ontime
636// and display it
637//
638TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
639{
640 TH1D collarea(area.GetHEnergy());
641 TH1D spectrum(excess);
642 TH1D weights;
643 hest.GetWeights(weights);
644
645 cout << "Effective On time: " << ontime << "s" << endl;
646
647 spectrum.SetDirectory(NULL);
648 spectrum.SetBit(kCanDelete);
649 spectrum.Scale(1./ontime);
650 spectrum.Divide(&collarea);
651 spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
652 spectrum.SetYTitle("N/sm^{2}");
653
654 TCanvas &c1 = fDisplay->AddTab("Spectrum");
655 c1.Divide(2,2);
656 c1.cd(1);
657 gPad->SetBorderMode(0);
658 gPad->SetLogx();
659 gPad->SetLogy();
660 gPad->SetGridx();
661 gPad->SetGridy();
662 collarea.DrawCopy();
663
664 c1.cd(2);
665 gPad->SetBorderMode(0);
666 gPad->SetLogx();
667 gPad->SetLogy();
668 gPad->SetGridx();
669 gPad->SetGridy();
670 TH1D *spec=(TH1D*)spectrum.DrawCopy();
671 FitSpectrum(*spec);
672
673 c1.cd(3);
674 gPad->SetBorderMode(0);
675 gPad->SetLogx();
676 gPad->SetLogy();
677 gPad->SetGridx();
678 gPad->SetGridy();
679 weights.DrawCopy();
680
681 //spectrum.Multiply(&weights);
682 spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
683 spectrum.SetBit(TH1::kNoStats);
684
685 for (int i=0; i<excess.GetNbinsX(); i++)
686 {
687 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
688 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
689
690 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
691 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
692 }
693
694 c1.cd(4);
695 gPad->SetBorderMode(0);
696 gPad->SetLogx();
697 gPad->SetLogy();
698 gPad->SetGridx();
699 gPad->SetGridy();
700 spectrum.SetMinimum(1e-12);
701 spectrum.SetXTitle("E [GeV]");
702 spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
703 spec = (TH1D*)spectrum.DrawCopy();
704
705 // Calculate Spectrum * E^2
706 for (int i=0; i<spectrum.GetNbinsX(); i++)
707 {
708 const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)*spectrum.GetBinLowEdge(i+2))*1e-3;
709
710 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
711 spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) *e*e);
712 }
713
714 // Display dN/dE*E^2 for conveinience
715 fDisplay->AddTab("Check");
716 gPad->SetLogx();
717 gPad->SetLogy();
718 gPad->SetGrid();
719
720 spectrum.SetName("FluxStd");
721 spectrum.SetMarkerStyle(kFullDotMedium);
722 spectrum.SetTitle("Differential flux times E^{2}");
723 spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
724 spectrum.SetDirectory(0);
725 spectrum.DrawCopy();
726
727 // Fit Spectrum
728 c1.cd(4);
729 return FitSpectrum(*spec);
730
731/*
732 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
733 f.SetParameter(0, -2.87);
734 f.SetParameter(1, 1.9e-6);
735 f.SetLineColor(kGreen);
736 spectrum.Fit(&f, "NIM", "", 100, 5000);
737 f.DrawCopy("same");
738
739 const Double_t p0 = f.GetParameter(0);
740 const Double_t p1 = f.GetParameter(1);
741
742 const Double_t e0 = f.GetParError(0);
743 const Double_t e1 = f.GetParError(1);
744
745 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
746 const Double_t exp = TMath::Power(10, np);
747
748 TString str;
749 str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
750 str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
751 str += "\\frac{ph}{TeVm^{2}s}";
752
753 TLatex tex;
754 tex.SetTextSize(0.045);
755 tex.SetBit(TLatex::kTextNDC);
756 tex.DrawLatex(0.45, 0.935, str);
757
758 str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
759 tex.DrawLatex(0.70, 0.83, str);
760
761 TArrayD res(2);
762 res[0] = f.GetParameter(0);
763 res[1] = f.GetParameter(1);
764
765 return res;
766 */
767/*
768 // Plot other spectra from Whipple
769 f.SetParameter(0, -2.45);
770 f.SetParameter(1, 3.3e-7);
771 f.SetRange(300, 8000);
772 f.SetLineColor(kBlack);
773 f.SetLineStyle(kDashed);
774 f.DrawCopy("same");
775
776 // Plot other spectra from Cangaroo
777 f.SetParameter(0, -2.53);
778 f.SetParameter(1, 2.0e-7);
779 f.SetRange(7000, 50000);
780 f.SetLineColor(kBlack);
781 f.SetLineStyle(kDashed);
782 f.DrawCopy("same");
783
784 // Plot other spectra from Robert
785 f.SetParameter(0, -2.59);
786 f.SetParameter(1, 2.58e-7);
787 f.SetRange(150, 1500);
788 f.SetLineColor(kBlack);
789 f.SetLineStyle(kDashed);
790 f.DrawCopy("same");
791
792 // Plot other spectra from HEGRA
793 f.SetParameter(0, -2.61);
794 f.SetParameter(1, 2.7e-7);
795 f.SetRange(1000, 20000);
796 f.SetLineColor(kBlack);
797 f.SetLineStyle(kDashed);
798 f.DrawCopy("same");
799 */
800}
801
802// --------------------------------------------------------------------------
803//
804// Scale some image parameter plots using the scale factor and plot them
805// together with the corresponding MC histograms.
806// Called from DisplaySize
807//
808Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
809{
810 TString same(name);
811 same += "Same";
812
813 TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
814 TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
815 if (!h1 || !h2)
816 return kFALSE;
817
818 TObject *obj = plist.FindObject(plot);
819 if (!obj)
820 {
821 *fLog << warn << plot << " not in parameter list... skipping." << endl;
822 return kFALSE;
823 }
824
825 TH1 *h3 = (TH1*)obj->FindObject(name);
826 if (!h3)
827 {
828 *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
829 return kFALSE;
830 }
831
832
833 const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
834 const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
835
836 gPad->SetBorderMode(0);
837 h2->SetLineColor(kBlack);
838 h3->SetLineColor(kBlue);
839 h2->Add(h1, -ascale);
840
841 //h2->Scale(1./ontime); //h2->Integral());
842 h3->Scale(scale); //h3->Integral());
843
844 h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
845
846 h2 = h2->DrawCopy();
847 h3 = h3->DrawCopy("same");
848
849 // Don't do this on the original object!
850 h2->SetStats(kFALSE);
851 h3->SetStats(kFALSE);
852
853 return kTRUE;
854}
855
856// --------------------------------------------------------------------------
857//
858// Take a lot of histograms and plot them together in one plot.
859// Calls PlotSame
860//
861Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
862{
863 *fLog << inf << "Reading from file: " << fPathIn << endl;
864
865 TFile file(fPathIn, "READ");
866 if (!file.IsOpen())
867 {
868 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
869 return kFALSE;
870 }
871
872 file.cd();
873 MStatusArray arr;
874 if (arr.Read()<=0)
875 {
876 *fLog << "MStatusDisplay not found in file... abort." << endl;
877 return kFALSE;
878 }
879
880 TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
881 if (!excess)
882 return kFALSE;
883
884 // ------------------- Plot excess versus size -------------------
885
886 TCanvas &c = fDisplay->AddTab("Excess");
887 c.Divide(3,2);
888 c.cd(1);
889 gPad->SetBorderMode(0);
890 gPad->SetLogx();
891 gPad->SetLogy();
892 gPad->SetGridx();
893 gPad->SetGridy();
894
895 excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
896 excess = excess->DrawCopy("E2");
897 // Don't do this on the original object!
898 excess->SetStats(kFALSE);
899 excess->SetMarkerStyle(kFullDotMedium);
900 excess->SetFillColor(kBlack);
901 excess->SetFillStyle(0);
902 excess->SetName("Excess ");
903 excess->SetDirectory(0);
904
905 TObject *o=0;
906 if ((o=plist.FindObject("ExcessMC")))
907 {
908 TH1 *histsel = (TH1F*)o->FindObject("");
909 if (histsel)
910 {
911 if (scale<0)
912 scale = excess->Integral()/histsel->Integral();
913
914 histsel->Scale(scale);
915 histsel->SetLineColor(kBlue);
916 histsel->SetBit(kCanDelete);
917 histsel = histsel->DrawCopy("E1 same");
918 // Don't do this on the original object!
919 histsel->SetStats(kFALSE);
920
921 fLog->Separator("Kolmogorov Test");
922 histsel->KolmogorovTest(excess, "DX");
923 fLog->Separator("Chi^2 Test");
924 const Double_t p = histsel->Chi2Test(excess, "P");
925
926 TLatex tex;
927 tex.SetBit(TLatex::kTextNDC);
928 tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
929 }
930 }
931
932 // -------------- Comparison of Image Parameters --------------
933 c.cd(2);
934 PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
935
936 c.cd(3);
937 PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
938
939 c.cd(4);
940 PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
941
942 c.cd(5);
943 PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
944
945 c.cd(6);
946 PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
947
948 return kTRUE;
949}
950
951Bool_t MJSpectrum::Process(const MDataSet &set)
952{
953 if (!set.IsValid())
954 {
955 *fLog << err << "ERROR - DataSet invalid!" << endl;
956 return kFALSE;
957 }
958
959 CheckEnv();
960
961 // --------------------------------------------------------------------------------
962
963 *fLog << inf;
964 fLog->Separator(GetDescriptor());
965 *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
966 *fLog << endl;
967
968 MBinning bins1("BinningAlpha");
969 MBinning bins2("BinningEnergyEst");
970 MBinning bins3("BinningTheta");
971 MBinning bins4("BinningFalseSource");
972 MBinning bins5("BinningWidth");
973 MBinning bins6("BinningLength");
974 MBinning bins7("BinningDist");
975 MBinning bins8("BinningMaxDist");
976 MBinning bins9("BinningM3Long");
977 MBinning bins0("BinningConc1");
978
979 MAlphaFitter fit;
980
981 MParList plist;
982 plist.AddToList(&bins1);
983 plist.AddToList(&bins3);
984 plist.AddToList(&bins4);
985 plist.AddToList(&bins5);
986 plist.AddToList(&bins6);
987 plist.AddToList(&bins7);
988 plist.AddToList(&bins8);
989 plist.AddToList(&bins9);
990 plist.AddToList(&bins0);
991 plist.AddToList(&fit);
992
993 TH1D temp1, size;
994 const Float_t ontime = ReadInput(plist, temp1, size);
995 if (ontime<0)
996 {
997 *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
998 return kFALSE;
999 }
1000
1001 plist.AddToList(&bins2);
1002
1003 MMcSpectrumWeight weight;
1004 if (!InitWeighting(set, weight))
1005 return kFALSE;
1006
1007 PrintSetup(fit);
1008 bins3.SetEdges(temp1, 'x');
1009
1010 TH1D temp2(temp1);
1011 if (!ReadOrigMCDistribution(set, temp2, weight))
1012 return kFALSE;
1013
1014 if (!GetThetaDistribution(temp1, temp2))
1015 return kFALSE;
1016
1017 if (!fNoThetaWeights)
1018 weight.SetWeightsZd(&temp1);
1019
1020 TH1D excess;
1021 if (!Refill(plist, excess))
1022 return kFALSE;
1023
1024 TH2D hist;
1025 MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
1026 if (fSimpleMode)
1027 {
1028 hist.UseCurrentStyle();
1029 MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
1030 if (!ReadOrigMCDistribution(set, hist, weight))
1031 return kFALSE;
1032
1033 if (!fRawMc)
1034 {
1035 for (int y=0; y<hist.GetNbinsY(); y++)
1036 for (int x=0; x<hist.GetNbinsX(); x++)
1037 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
1038 //hist.SetEntries(hist.Integral());
1039 }
1040 }
1041 else
1042 {
1043 weight.SetNameMcEvt("MMcEvtBasic");
1044 if (!IntermediateLoop(plist, mh1, temp1, set, weight))
1045 return kFALSE;
1046 weight.SetNameMcEvt();
1047 }
1048
1049 DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1050
1051 // ------------------------- Final loop --------------------------
1052
1053 *fLog << endl;
1054 fLog->Separator("Calculate efficiencies");
1055 *fLog << endl;
1056
1057 MTaskList tlist2;
1058 plist.Replace(&tlist2);
1059
1060 MReadMarsFile read("Events");
1061 read.DisableAutoScheme();
1062 set.AddFilesOn(read);
1063
1064 // Selector to get correct (final) theta-distribution
1065 temp1.SetXTitle("MPointingPos.fZd");
1066
1067 MH3 mh3(temp1);
1068
1069 MFEventSelector2 sel2(mh3);
1070 sel2.SetHistIsProbability();
1071
1072 MContinue contsel(&sel2);
1073 contsel.SetInverted();
1074
1075 // Get correct source position
1076 //MSrcPosCalc calc;
1077
1078 // Calculate corresponding Hillas parameters
1079 /*
1080 MHillasCalc hcalc1;
1081 MHillasCalc hcalc2("MHillasCalcAnti");
1082 hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1083 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1084 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1085 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1086 */
1087
1088 // Fill collection area and energy estimator (unfolding)
1089 // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1090 MHCollectionArea area0("TriggerArea");
1091 MHCollectionArea area1;
1092 area0.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1093 area1.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1094
1095 MHEnergyEst hest;
1096
1097 MFillH fill30(&area0, "", "FillTriggerArea");
1098 MFillH fill31(&area1, "", "FillCollectionArea");
1099 MFillH fill4(&hest, "", "FillEnergyEst");
1100 MFillH fill5("MHThreshold", "", "FillThreshold");
1101 fill30.SetWeight();
1102 fill31.SetWeight();
1103 fill4.SetWeight();
1104 fill5.SetWeight();
1105 fill30.SetNameTab("TrigArea");
1106 fill31.SetNameTab("ColArea");
1107 fill4.SetNameTab("E-Est");
1108 fill5.SetNameTab("Threshold");
1109
1110 MH3 hsize("MHillas.fSize");
1111 hsize.SetName("ExcessMC");
1112 hsize.Sumw2();
1113
1114 MBinning bins(size, "BinningExcessMC");
1115 plist.AddToList(&hsize);
1116 plist.AddToList(&bins);
1117
1118 MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1119 MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1120 MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1121 MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1122 MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1123 MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1124 MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1125 MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1126 fill1a.SetNameTab("PreCut");
1127 fill2a.SetNameTab("PostCut");
1128 fill3a.SetNameTab("VsSize");
1129 fill4a.SetNameTab("HilExt");
1130 fill5a.SetNameTab("HilSrc");
1131 fill6a.SetNameTab("ImgPar");
1132 fill7a.SetNameTab("NewPar");
1133 fill8a.SetBit(MFillH::kDoNotDisplay);
1134 fill1a.SetWeight();
1135 fill2a.SetWeight();
1136 fill3a.SetWeight();
1137 fill4a.SetWeight();
1138 fill5a.SetWeight();
1139 fill6a.SetWeight();
1140 fill7a.SetWeight();
1141 fill8a.SetWeight();
1142
1143 MEnergyEstimate est;
1144 MTaskEnv taskenv1("EstimateEnergy");
1145 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1146
1147 tlist2.AddToList(&read);
1148 if (!fRawMc && fNoThetaWeights)
1149 tlist2.AddToList(&contsel);
1150 //tlist2.AddToList(&calc);
1151 //tlist2.AddToList(&hcalc1);
1152 //tlist2.AddToList(&hcalc2);
1153 tlist2.AddToList(&weight);
1154 tlist2.AddToList(&fill1a);
1155 tlist2.AddToList(&fill30);
1156 tlist2.AddToList(fCut0);
1157 tlist2.AddToList(fCut1);
1158 tlist2.AddToList(fCut2);
1159 tlist2.AddToList(fCut3);
1160 tlist2.AddToList(&taskenv1);
1161 tlist2.AddToList(&fill31);
1162 tlist2.AddToList(&fill4);
1163 tlist2.AddToList(&fill5);
1164 tlist2.AddToList(&fill2a);
1165 tlist2.AddToList(&fill3a);
1166 tlist2.AddToList(&fill4a);
1167 tlist2.AddToList(&fill5a);
1168 tlist2.AddToList(&fill6a);
1169 tlist2.AddToList(&fill7a);
1170 tlist2.AddToList(&fill8a);
1171 //tlist2.AddToList(&fill9a);
1172
1173 // by setting it here it is distributed to all consecutive tasks
1174 tlist2.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
1175
1176 MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
1177 loop2.SetParList(&plist);
1178 loop2.SetDisplay(fDisplay);
1179 loop2.SetLogStream(fLog);
1180
1181 if (!SetupEnv(loop2))
1182 return kFALSE;
1183
1184 if (!loop2.Eventloop(fMaxEvents))
1185 {
1186 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1187 return kFALSE;
1188 }
1189
1190 if (!loop2.GetDisplay())
1191 {
1192 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1193 return kFALSE;
1194 }
1195
1196 gLog.Separator("Energy Estimator");
1197 if (plist.FindObject("EstimateEnergy"))
1198 plist.FindObject("EstimateEnergy")->Print();
1199
1200 gLog.Separator("Spectrum");
1201
1202 // -------------------------- Spectrum ----------------------------
1203
1204 // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1205 TArrayD res(DisplaySpectrum(area1, excess, hest, ontime));
1206
1207 // Spectrum fitted (convert res[1] from TeV to GeV)
1208 TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
1209
1210 // Number of events this spectrum would produce per s and m^2
1211 Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1212
1213 // scale with effective collection area to get the event rate (N/s)
1214 // scale with the effective observation time to absolute observed events
1215 n *= area1.GetCollectionAreaAbs()*ontime; // N
1216
1217 // Now calculate the scale factor from the number of events
1218 // produced and the number of events which should have been
1219 // observed with our telescope in the time ontime
1220 const Double_t scale = n/area1.GetEntries();
1221
1222 // Print normalization constant
1223 cout << "MC normalization factor: " << scale << endl;
1224
1225 // Overlay normalized plots
1226 DisplaySize(plist, scale);
1227
1228 // check if output should be written
1229 if (!fPathOut.IsNull())
1230 fDisplay->SaveAsRoot(fPathOut);
1231
1232 return kTRUE;
1233}
Note: See TracBrowser for help on using the repository browser.