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

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