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

Last change on this file since 8722 was 8719, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 49.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-2007
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 "MDirIter.h"
50
51#include "MStatusArray.h"
52#include "MStatusDisplay.h"
53
54// Container
55#include "MH3.h"
56#include "MHn.h"
57#include "MBinning.h"
58#include "MDataSet.h"
59#include "MMcCorsikaRunHeader.h"
60
61// Spectrum
62#include "../mhflux/MAlphaFitter.h"
63#include "../mhflux/MHAlpha.h"
64#include "../mhflux/MHCollectionArea.h"
65#include "../mhflux/MHEnergyEst.h"
66#include "../mhflux/MMcSpectrumWeight.h"
67
68// Eventloop
69#include "MEvtLoop.h"
70#include "MTaskList.h"
71#include "MParList.h"
72
73// Tasks/Filter
74#include "MReadMarsFile.h"
75#include "MEnergyEstimate.h"
76#include "MTaskEnv.h"
77#include "MFillH.h"
78#include "MFDataPhrase.h"
79#include "MFDataMember.h"
80#include "MContinue.h"
81
82ClassImp(MJSpectrum);
83
84using namespace std;
85
86MJSpectrum::MJSpectrum(const char *name, const char *title)
87 : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0),
88 fEstimateEnergy(0), fCalcHadronness(0), fCalcDisp(0), fForceTheta(kFALSE)
89{
90 fName = name ? name : "MJSpectrum";
91 fTitle = title ? title : "Standard program to calculate spectrum";
92
93 // Make sure that fDisplay is maintained properly
94 // (i.e. removed via RecursiveRemove if closed)
95 gROOT->GetListOfCleanups()->Add(this);
96}
97
98// --------------------------------------------------------------------------
99//
100// Delete all the fCut* data members and fCalcHadronness/fEstimateEnergy
101//
102MJSpectrum::~MJSpectrum()
103{
104 if (fCutQ)
105 delete fCutQ;
106 if (fCut0)
107 delete fCut0;
108 if (fCut1)
109 delete fCut1;
110 if (fCut2)
111 delete fCut2;
112 if (fCut3)
113 delete fCut3;
114 if (fCutS)
115 delete fCutS;
116 if (fEstimateEnergy)
117 delete fEstimateEnergy;
118 if (fCalcHadronness)
119 delete fCalcHadronness;
120 if (fCalcDisp)
121 delete fCalcDisp;
122}
123
124// --------------------------------------------------------------------------
125//
126// Setup a task estimating the energy. The given task is cloned.
127//
128void MJSpectrum::SetEnergyEstimator(const MTask *task)
129{
130 if (fEstimateEnergy)
131 delete fEstimateEnergy;
132 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
133}
134
135// --------------------------------------------------------------------------
136//
137// Read a MTask named name into task from the open file. If this task is
138// required mustexist can be set. Depending on success kTRUE or kFALSE is
139// returned. If the task is a MContinue SetAllowEmpty is called.
140// The name of the task is set to name.
141//
142Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const
143{
144 // If a task is set delete task
145 if (task)
146 {
147 delete task;
148 task = 0;
149 }
150
151 // Try to get task from file
152 task = (MTask*)gFile->Get(name);
153 if (!task)
154 {
155 if (!mustexist)
156 return kTRUE;
157 *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
158 return kFALSE;
159 }
160
161 // Check if task inherits from MTask
162 if (!task->InheritsFrom(MTask::Class()))
163 {
164 *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
165 delete task;
166 return kFALSE;
167 }
168
169 // Set name of task
170 task->SetName(name);
171
172 // SetAllowEmpty for MContinue tasks
173 if (dynamic_cast<MContinue*>(task))
174 dynamic_cast<MContinue*>(task)->SetAllowEmpty();
175
176 return kTRUE;
177}
178
179// --------------------------------------------------------------------------
180//
181// Print setup used for the MC processing, namely MAlphaFitter ans all fCut*.
182//
183void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
184{
185 fLog->Separator("Alpha Fitter");
186 *fLog << all;
187 fit.Print("result");
188
189 fLog->Separator("Used Cuts");
190 fCutQ->Print();
191 fCut0->Print();
192 fCut1->Print();
193 fCutS->Print();
194 fCut2->Print();
195 fCut3->Print();
196}
197
198// --------------------------------------------------------------------------
199//
200// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
201// the dataset.
202// The simulated energy range and spectral slope is initialized from
203// there.
204// In the following eventloops the forced check in MMcSpectrumWeight
205// ensures, that the spectral slope and energy range doesn't change.
206//
207Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
208{
209 fLog->Separator("Initialize energy weighting");
210
211 // Read Resources
212 if (!CheckEnv(w))
213 {
214 *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
215 return kFALSE;
216 }
217
218 w.Print("new");
219
220 return kTRUE;
221}
222
223Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
224{
225 *fLog << inf << "Reading from file: " << fPathIn << endl;
226
227 TFile file(fPathIn, "READ");
228 if (!file.IsOpen())
229 {
230 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
231 return -1;
232 }
233
234 MStatusArray arr;
235 if (arr.Read()<=0)
236 {
237 *fLog << "MStatusDisplay not found in file... abort." << endl;
238 return -1;
239 }
240
241 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
242 TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
243 if (!vstime || !size)
244 return -1;
245
246 vstime->Copy(h1);
247 size->Copy(h2);
248 h1.SetDirectory(0);
249 h2.SetDirectory(0);
250
251 if (fDisplay)
252 arr.DisplayIn(*fDisplay, "Hist");
253
254 if (!ReadTask(fCutQ, "CutQ"))
255 return -1;
256 if (!ReadTask(fCut0, "Cut0"))
257 return -1;
258 if (!ReadTask(fCut1, "Cut1"))
259 return -1;
260 if (!ReadTask(fCut2, "Cut2"))
261 return -1;
262 if (!ReadTask(fCut3, "Cut3"))
263 return -1;
264 if (!ReadTask(fCalcHadronness, "CalcHadronness", kFALSE))
265 return -1;
266 if (!ReadTask(fCalcDisp, "CalcDisp", kFALSE))
267 return -1;
268
269 TObjArray arrread;
270
271 TIter Next(plist);
272 TObject *o=0;
273 while ((o=Next()))
274 if (o->InheritsFrom(MBinning::Class()))
275 arrread.Add(o);
276
277 arrread.Add(plist.FindObject("MAlphaFitter"));
278
279 if (!ReadContainer(arrread))
280 return -1;
281
282 return vstime->Integral();
283}
284
285// --------------------------------------------------------------------------
286//
287// return maximum value of MMcRunHeader.fImpactMax stored in the RunHeaders
288// of the files from the MC dataset
289//
290Bool_t MJSpectrum::AnalyzeMC(const MDataSet &set, Float_t &impactmax, Float_t &emin/*, Float_t emax*/) const
291{
292 if (fDisplay)
293 fDisplay->SetStatusLine1("Analyzing Monte Carlo headers...");
294
295 // ----- Create chain -----
296 *fLog << inf << "Getting files... " << flush;
297 TChain chain("RunHeaders");
298 if (!set.AddFilesOn(chain))
299 return kFALSE;
300 *fLog << "done. " << endl;
301
302 *fLog << all;
303 *fLog << "Searching for maximum impact... " << flush;
304 impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax");
305 *fLog << "found " << impactmax/100 << "m" << endl;
306
307 *fLog << "Searching for minimum lower energy limit... " << flush;
308 emin = chain.GetMinimum("MMcCorsikaRunHeader.fELowLim");
309 *fLog << "found " << emin << "GeV" << endl;
310
311 *fLog << "Searching for maximum lower energy limit... " << flush;
312 const Float_t emin2 = chain.GetMaximum("MMcCorsikaRunHeader.fELowLim");
313 *fLog << "found " << emin2 << "GeV" << endl;
314
315 if (emin2>emin)
316 {
317 *fLog << warn;
318 *fLog << "WARNING - You are using files with different lower limits for the simulated" << endl;
319 *fLog << " energy, i.e. that the collection area and your correction" << endl;
320 *fLog << " coefficients between " << emin << "GeV and ";
321 *fLog << emin2 << "GeV might be wrong." << endl;
322 }
323
324/*
325 *fLog << "Getting maximum energy... " << flush;
326 emax = chain.GetMaximum("MMcCorsikaRunHeader.fEUppLim");
327 *fLog << "found " << emax << "GeV" << endl;
328 */
329 return kTRUE;
330}
331
332Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
333{
334 // Some debug output
335 fLog->Separator("Compiling original MC distribution");
336
337 Float_t impactmax=0, Emin=0;
338 if (!AnalyzeMC(set, impactmax, Emin))
339 return kFALSE;
340
341 *fLog << inf << "Getting files... " << flush;
342 MDirIter iter;
343 if (!set.AddFilesOn(iter))
344 return kFALSE;
345 *fLog << "done. " << endl;
346
347 const Int_t tot = iter.GetNumEntries();
348
349 // Prepare histogram
350 h.Reset();
351 h.Sumw2();
352 if (h.InheritsFrom(TH2::Class()))
353 {
354 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
355 h.SetXTitle("\\Theta [\\circ]");
356 h.SetYTitle("E [GeV]");
357 h.SetZTitle(Form("Counts normalized to r=%.0fm", impactmax/100));
358 }
359 else
360 {
361 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
362 h.SetXTitle("\\Theta [\\circ]");
363 h.SetYTitle(Form("Counts normalized to r=%.0fm", impactmax/100));
364 }
365 h.SetDirectory(0);
366
367 const TString cont = h.InheritsFrom(TH2::Class()) ?
368 "MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC" :
369 "MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC";
370
371 if (fDisplay)
372 fDisplay->SetStatusLine1("Reading OriginalMC distribution...");
373
374 TH1 *hfill = (TH1*)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC");
375 hfill->SetDirectory(0);
376
377 *fLog << all << "Compiling simulated source spectrum... please stand by, this may take a while." << endl;
378
379 Double_t evts = 0;
380 Int_t num = 0;
381
382 // Reading this with a eventloop is five times slower :(
383 TString fname;
384 while (fDisplay)
385 {
386 if (fDisplay)
387 fDisplay->SetProgressBarPosition(Float_t(num++)/tot);
388
389 // Get next filename
390 fname = iter.Next();
391 if (fname.IsNull())
392 break;
393
394 if (fDisplay)
395 fDisplay->SetStatusLine2(fname);
396
397 // open file
398 TFile file(fname);
399 if (file.IsZombie())
400 {
401 *fLog << err << "ERROR - Couldn't open file " << fname << endl;
402 return kFALSE;
403 }
404
405 // Get tree OriginalMC
406 TTree *t = dynamic_cast<TTree*>(file.Get("OriginalMC"));
407 if (!t)
408 {
409 *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl;
410 return kFALSE;
411 }
412 if (t->GetEntries()==0)
413 {
414 *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl;
415 continue;
416 }
417
418 // Get tree RunHeaders
419 TTree *rh = dynamic_cast<TTree*>(file.Get("RunHeaders"));
420 if (!rh)
421 {
422 *fLog << err << "ERROR - File " << fname << " doesn't contain tree RunHeaders." << endl;
423 return kFALSE;
424 }
425 if (rh->GetEntries()!=1)
426 {
427 *fLog << err << "ERROR - RunHeaders of " << fname << " doesn't contain exactly one entry." << endl;
428 return kFALSE;
429 }
430
431 // Get corsika run header
432 MMcCorsikaRunHeader *head=0;
433 rh->SetBranchAddress("MMcCorsikaRunHeader.", &head);
434 rh->GetEntry(0);
435 if (!head)
436 {
437 *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from " << fname << "." << endl;
438 return kFALSE;
439 }
440
441 // Get the maximum impact parameter of this file. Due to different
442 // production areas an additional scale-factor is applied.
443 //
444 // Because it is assumed that the efficiency outside the production
445 // area is nearly zero no additional weight has to be applied to the
446 // events after cuts.
447 const Double_t impact = rh->GetMaximum("MMcRunHeader.fImpactMax");
448 const Double_t scale = impactmax/impact;
449
450 // Propagate the run header to MMcSpectrumWeight
451 if (!weight.Set(*head))
452 {
453 *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
454 return kFALSE;
455 }
456
457 // Get the corresponding rule for the weights
458 const TString weights(weight.GetFormulaWeights("MMcEvtBasic"));
459
460 // No we found everything... go on reading contents
461 *fLog << inf << "Reading OriginalMC of " << fname << endl;
462
463 // Fill histogram from tree
464 hfill->SetDirectory(&file);
465 if (t->Draw(cont, weights, "goff")<0)
466 {
467 *fLog << err << "ERROR - Reading OriginalMC failed." << endl;
468 return kFALSE;
469 }
470 hfill->SetDirectory(0);
471
472 evts += hfill->GetEntries();
473
474 // ----- Complete incomplete energy ranges -----
475 // ----- and apply production area weights -----
476 weight.CompleteEnergySpectrum(*hfill, Emin);
477
478 hfill->Scale(scale*scale);
479
480 // Add weighted data from file to final histogram
481 h.Add(hfill);
482 }
483 delete hfill;
484
485 *fLog << all << "Total number of events from files in Monte Carlo dataset: " << evts << endl;
486 if (fDisplay)
487 fDisplay->SetStatusLine2("done.");
488
489 if (!fDisplay || h.GetEntries()>0)
490 return kTRUE;
491
492 *fLog << err << "ERROR - Histogram with distribution from OriginalMC empty..." << endl;
493 return kFALSE;
494}
495
496void MJSpectrum::GetThetaDistribution(TH1D &temp1, TH2D &hist2) const
497{
498 TH1D &temp2 = *hist2.ProjectionX();
499
500 // Display some stuff
501 if (fDisplay)
502 {
503 TCanvas &c = fDisplay->AddTab("ZdDist");
504 c.Divide(2,2);
505
506 // On-Time vs. Theta
507 c.cd(1);
508 gPad->SetBorderMode(0);
509 gPad->SetGridx();
510 gPad->SetGridy();
511 temp1.DrawCopy();
512
513 // Number of MC events (produced) vs Theta
514 c.cd(2);
515 gPad->SetBorderMode(0);
516 gPad->SetGridx();
517 gPad->SetGridy();
518 temp2.SetName("NVsTheta");
519 temp2.DrawCopy("hist");
520
521 c.cd(4);
522 gPad->SetBorderMode(0);
523 gPad->SetGridx();
524 gPad->SetGridy();
525
526 c.cd(3);
527 gPad->SetBorderMode(0);
528 gPad->SetGridx();
529 gPad->SetGridy();
530 }
531
532 // Calculate the Probability
533 temp1.Divide(&temp2);
534 temp1.Scale(1./temp1.Integral());
535
536 // Some cosmetics: Name, Axis, etc.
537 temp1.SetName("ProbVsTheta");
538 temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
539 temp1.SetYTitle("Probability");
540 if (fDisplay)
541 temp1.DrawCopy("hist");
542
543 delete &temp2;
544}
545
546// --------------------------------------------------------------------------
547//
548// Display the final theta distribution.
549//
550Bool_t MJSpectrum::DisplayResult(const TH2D &h2) const
551{
552 if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
553 {
554 *fLog << err << "ERROR - Display or tab ZdDist vanished... abort." << endl;
555 return kFALSE;
556 }
557
558 TH1D &proj = *h2.ProjectionX();
559 proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
560 proj.SetXTitle("\\Theta [\\circ]");
561 proj.SetYTitle("Counts");
562 proj.SetLineColor(kBlue);
563 proj.SetDirectory(0);
564 proj.SetBit(kCanDelete);
565
566 TVirtualPad *pad = gPad;
567
568 pad->cd(4);
569 proj.DrawCopy();
570
571 pad->cd(1);
572 TH1D *theta = (TH1D*)gPad->FindObject("Theta");
573 if (!theta)
574 {
575 *fLog << err << "ERROR - Theta-Histogram vanished... cannot proceed." << endl;
576 return kFALSE;
577 }
578
579 // Check whether histogram is empty
580 if (proj.GetMaximum()==0)
581 {
582 *fLog << err;
583 *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap" << endl;
584 *fLog << " with the Zenith Angle distribution of your observation." << endl;
585 *fLog << " Maybe the energy binning is undefined or wrong (from ";
586 *fLog << h2.GetYaxis()->GetXmin() << "GeV to " << h2.GetYaxis()->GetXmax() << "GeV)" << endl;
587 theta->SetLineColor(kRed);
588 return kFALSE;;
589 }
590
591 // scale histogram and set new maximum for display
592 proj.Scale(theta->GetMaximum()/proj.GetMaximum());
593 theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
594
595 // draw project
596 proj.Draw("same");
597
598 // Compare both histograms
599 *fLog << inf << "Comparing theta distributions for data and MCs." << endl;
600
601 const Double_t prob = proj.Chi2Test(theta, "P");
602 if (prob==1)
603 return kTRUE;
604
605 if (prob>0.99)
606 {
607 *fLog << inf;
608 *fLog << "The Zenith Angle distribution of your Monte Carlos fits well" << endl;
609 *fLog << "with the Zenith Angle distribution of your observation." << endl;
610 *fLog << "The probability for identical Theta distributions is " << prob << endl;
611 return kTRUE;
612 }
613
614 if (prob<0.01)
615 {
616 *fLog << err;
617 *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit" << endl;
618 *fLog << " with the Zenith Angle distribution of your observation." << endl;
619 *fLog << " The probability for identical Theta distributions is " << prob << endl;
620 if (!fForceTheta)
621 *fLog << " To force processing use --force-theta (with care!)" << endl;
622 theta->SetLineColor(kRed);
623 return fForceTheta;
624 }
625
626 *fLog << warn;
627 *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well" << endl;
628 *fLog << " with the Zenith Angle distribution of your observation." << endl;
629 *fLog << " The probability for identical Theta distributions is " << prob << endl;
630 return kTRUE;
631}
632
633// --------------------------------------------------------------------------
634//
635// Fills the excess histogram (vs E-est) from the events stored in the
636// ganymed result file and therefor estimates the energy.
637//
638// The resulting histogram excess-vs-energy ist copied into h2.
639//
640Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2)/*const*/
641{
642 // Try to find the class used to determin the signal!
643 TString cls("MHAlpha");
644 if (fDisplay)
645 {
646 TCanvas *c = fDisplay->GetCanvas("Hist");
647 if (c)
648 {
649 TIter Next(c->GetListOfPrimitives());
650 TObject *obj=0;
651 while ((obj=Next()))
652 if (obj->InheritsFrom(MHAlpha::Class()))
653 break;
654 if (obj)
655 cls = obj->ClassName();
656 }
657 }
658
659 // Now fill the histogram
660 *fLog << endl;
661 fLog->Separator("Refill Excess");
662 *fLog << endl;
663
664 MTaskList tlist;
665 plist.AddToList(&tlist);
666
667 MReadTree read("Events");
668 read.DisableAutoScheme();
669 read.AddFile(fPathIn);
670/*
671 MTaskEnv taskenv0("CalcDisp");
672 taskenv0.SetDefault(fCalcDisp);
673
674 MTaskEnv taskenv1("CalcHadronness");
675 taskenv1.SetDefault(fCalcHadronness);
676 */
677 MEnergyEstimate est;
678 MTaskEnv taskenv2("EstimateEnergy");
679 taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
680
681 MContinue *cont = new MContinue("", "CutS");
682 cont->SetAllowEmpty();
683
684 if (fCutS)
685 delete fCutS;
686 fCutS = cont;
687
688 // FIXME: Create HistE and HistEOff to be able to modify it from
689 // the resource file.
690
691 MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
692 MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
693
694 MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
695 MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
696
697 fill1.SetFilter(&f0);
698 fill2.SetFilter(&f1);
699
700 tlist.AddToList(&read);
701 //tlist.AddToList(&taskenv0); // not necessary, stored in file!
702 //tlist.AddToList(&taskenv1); // not necessary, stored in file!
703 tlist.AddToList(fCutS);
704 tlist.AddToList(&taskenv2);
705 tlist.AddToList(&f0);
706 tlist.AddToList(&f1);
707 tlist.AddToList(&fill1);
708 tlist.AddToList(&fill2);
709
710 // by setting it here it is distributed to all consecutive tasks
711 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
712
713 MEvtLoop loop("RefillExcess"); // ***** fName *****
714 loop.SetParList(&plist);
715 loop.SetDisplay(fDisplay);
716 loop.SetLogStream(fLog);
717
718 if (!SetupEnv(loop))
719 return kFALSE;
720
721 if (!loop.Eventloop())
722 {
723 *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
724 return kFALSE;
725 }
726
727 if (!loop.GetDisplay())
728 {
729 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
730 return kFALSE;
731 }
732
733 const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
734 if (!halpha)
735 {
736 *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
737 return kFALSE;
738 }
739
740 halpha->GetHEnergy().Copy(h2);
741 h2.SetDirectory(0);
742
743 return kTRUE;
744}
745
746/*
747Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
748{
749 MTaskList tlist1;
750 plist.Replace(&tlist1);
751
752 MReadMarsFile readmc("OriginalMC");
753 //readmc.DisableAutoScheme();
754 if (!set.AddFilesOn(readmc))
755 return kFALSE;
756
757 readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
758 readmc.EnableBranch("MMcEvtBasic.fEnergy");
759
760 mh1.SetLogy();
761 mh1.SetLogz();
762 mh1.SetName("ThetaE");
763
764 MFillH fill0(&mh1);
765 //fill0.SetDrawOption("projx only");
766
767 MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
768 MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
769 if (bins2 && bins3)
770 {
771 bins2->SetName("BinningThetaEY");
772 bins3->SetName("BinningThetaEX");
773 }
774 tlist1.AddToList(&readmc);
775 tlist1.AddToList(&weight);
776
777 temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
778 MH3 mh3mc(temp1);
779
780 MFEventSelector2 sel1(mh3mc);
781 sel1.SetHistIsProbability();
782
783 fill0.SetFilter(&sel1);
784
785 //if (!fRawMc)
786 tlist1.AddToList(&sel1);
787 tlist1.AddToList(&fill0);
788
789 // by setting it here it is distributed to all consecutive tasks
790 tlist1.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
791
792 MEvtLoop loop1("IntermediateLoop"); // ***** fName *****
793 loop1.SetParList(&plist);
794 loop1.SetLogStream(fLog);
795 loop1.SetDisplay(fDisplay);
796
797 if (!SetupEnv(loop1))
798 return kFALSE;
799
800 if (!loop1.Eventloop(fMaxEvents))
801 {
802 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
803 return kFALSE;
804 }
805
806 if (!loop1.GetDisplay())
807 {
808 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
809 return kFALSE;
810 }
811
812 if (bins2 && bins3)
813 {
814 bins2->SetName("BinningEnergyEst");
815 bins3->SetName("BinningTheta");
816 }
817
818 return kTRUE;
819}
820*/
821
822TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
823{
824 const Double_t p0 = f.GetParameter(0);
825 const Double_t p1 = f.GetParameter(1);
826
827 const Double_t e0 = f.GetParError(0);
828 const Double_t e1 = f.GetParError(1);
829
830 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
831 const Double_t exp = TMath::Power(10, np);
832
833 const Float_t fe0 = TMath::Log10(e0);
834 const Float_t fe1 = TMath::Log10(e1);
835
836 // calc number of (gueltige ziffern) digits to be displayed
837 const char *f0 = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0";
838 const char *f1 = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0";
839
840 TString str;
841 switch (type)
842 {
843 case 0:
844 {
845 const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1);
846 const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
847
848 str = Form(fmt0.Data(), p1/exp, e1/exp, np);
849 str += Form(fmt1.Data(), p0, e0);
850 str += "\\frac{ph}{TeVm^{2}s}";
851 }
852 break;
853 case 1:
854 str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
855 break;
856 case 2:
857 str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF()));
858 break;
859 }
860 return str;
861}
862
863TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
864{
865 Axis_t lo, hi;
866 MH::GetRangeUser(spectrum, lo, hi);
867
868 TF1 f("f", "[1]*(x/500)^[0]", lo, hi);
869 f.SetParameter(0, -3.0);
870 f.SetParameter(1, spectrum.GetMaximum());
871 f.SetLineColor(kBlue);
872 f.SetLineWidth(2);
873 spectrum.Fit(&f, "NIR"); // M skipped
874 f.DrawCopy("same");
875
876 TString str = FormString(f);
877
878 TLatex tex;
879 tex.SetTextSize(0.045);
880 tex.SetBit(TLatex::kTextNDC);
881 tex.SetTextAlign(31);
882 tex.DrawLatex(0.89, 0.935, str);
883
884 str = FormString(f, 1);
885 tex.DrawLatex(0.89, 0.83, str);
886
887 str = FormString(f, 2);
888 tex.DrawLatex(0.89, 0.735, str);
889
890 TArrayD res(2);
891 res[0] = f.GetParameter(0);
892 res[1] = f.GetParameter(1);
893 return res;
894}
895
896// --------------------------------------------------------------------------
897//
898// Calculate the final spectrum from:
899// - collection area
900// - excess
901// - correction coefficients
902// - ontime
903// and display it
904//
905TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
906{
907 // Create copies of the histograms
908 TH1D collarea(area.GetHEnergy());
909 TH1D spectrum(excess);
910 TH1D weights;
911
912 // Get spill-over corrections from energy estimation
913 hest.GetWeights(weights);
914
915 // Print effective on-time
916 cout << "Effective On time: " << ontime << "s" << endl;
917
918 // Setup spectrum plot
919 spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
920 spectrum.SetYTitle("N/sm^{2}");
921 spectrum.SetDirectory(NULL);
922 spectrum.SetBit(kCanDelete);
923
924 // Divide by collection are and on-time
925 spectrum.Scale(1./ontime);
926 spectrum.Divide(&collarea);
927
928 // Draw spectrum before applying spill-over corrections
929 TCanvas &c1 = fDisplay->AddTab("Spectrum");
930 c1.Divide(2,2);
931 c1.cd(1);
932 gPad->SetBorderMode(0);
933 gPad->SetLogx();
934 gPad->SetLogy();
935 gPad->SetGridx();
936 gPad->SetGridy();
937 collarea.DrawCopy();
938
939 c1.cd(2);
940 gPad->SetBorderMode(0);
941 gPad->SetLogx();
942 gPad->SetLogy();
943 gPad->SetGridx();
944 gPad->SetGridy();
945 TH1D *spec=(TH1D*)spectrum.DrawCopy();
946 FitSpectrum(*spec);
947
948 c1.cd(3);
949 gPad->SetBorderMode(0);
950 gPad->SetLogx();
951 gPad->SetLogy();
952 gPad->SetGridx();
953 gPad->SetGridy();
954 weights.DrawCopy();
955
956 // Apply spill-over correction (done't take the errors into account)
957 // They are supposed to be identical with the errors of the
958 // collection area and cancel out.
959 //spectrum.Multiply(&weights);
960 spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
961 spectrum.SetBit(TH1::kNoStats);
962
963 for (int i=0; i<excess.GetNbinsX(); i++)
964 {
965 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
966 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
967
968 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
969 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
970 }
971
972 c1.cd(4);
973 gPad->SetBorderMode(0);
974 gPad->SetLogx();
975 gPad->SetLogy();
976 gPad->SetGridx();
977 gPad->SetGridy();
978 spectrum.SetMinimum(1e-12);
979 spectrum.SetXTitle("E [GeV]");
980 spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
981 spec = (TH1D*)spectrum.DrawCopy();
982
983 TF1 fc("Crab", "7.0e-6*(x/300)^(-2.31-0.26*log10(x/300))", 100, 6000);
984 fc.SetLineStyle(7);
985 fc.SetLineWidth(2);
986 fc.SetLineColor(14);
987 fc.DrawCopy("same");
988
989 TF1 fa("PG1553", "1.8e-6*(x/200)^-4.21", 90, 600);
990 static_cast<const TAttLine&>(fc).Copy(fa);
991 fa.DrawCopy("same");
992
993 // Display dN/dE*E^2 for conveinience
994 fDisplay->AddTab("Check");
995 gPad->SetLogx();
996 gPad->SetLogy();
997 gPad->SetGrid();
998
999 // Calculate Spectrum * E^2
1000 for (int i=0; i<spectrum.GetNbinsX(); i++)
1001 {
1002 const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)*spectrum.GetBinLowEdge(i+2))*1e-3;
1003
1004 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
1005 spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) *e*e);
1006 }
1007
1008 spectrum.SetName("FluxStd");
1009 spectrum.SetMarkerStyle(kFullDotMedium);
1010 spectrum.SetTitle("Differential flux times E^{2}");
1011 spectrum.SetYTitle("E^{2}·dN/dE [N·TeV/sm^{2}]");
1012 spectrum.SetDirectory(0);
1013 spectrum.DrawCopy();
1014
1015 TF1 fc2("Crab*E^2", "x^2*Crab*1e-6", 100, 6000);
1016 static_cast<const TAttLine&>(fc).Copy(fc2);
1017 fc2.DrawCopy("same");
1018
1019 TF1 fa2("PG1553*E^2", "x^2*PG1553*1e-6", 100, 6000);
1020 static_cast<const TAttLine&>(fc).Copy(fa2);
1021 fa2.DrawCopy("same");
1022
1023 // Fit Spectrum
1024 c1.cd(4);
1025 return FitSpectrum(*spec);
1026
1027/*
1028 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
1029 f.SetParameter(0, -2.87);
1030 f.SetParameter(1, 1.9e-6);
1031 f.SetLineColor(kGreen);
1032 spectrum.Fit(&f, "NIM", "", 100, 5000);
1033 f.DrawCopy("same");
1034
1035 const Double_t p0 = f.GetParameter(0);
1036 const Double_t p1 = f.GetParameter(1);
1037
1038 const Double_t e0 = f.GetParError(0);
1039 const Double_t e1 = f.GetParError(1);
1040
1041 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
1042 const Double_t exp = TMath::Power(10, np);
1043
1044 TString str;
1045 str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
1046 str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
1047 str += "\\frac{ph}{TeVm^{2}s}";
1048
1049 TLatex tex;
1050 tex.SetTextSize(0.045);
1051 tex.SetBit(TLatex::kTextNDC);
1052 tex.DrawLatex(0.45, 0.935, str);
1053
1054 str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
1055 tex.DrawLatex(0.70, 0.83, str);
1056
1057 TArrayD res(2);
1058 res[0] = f.GetParameter(0);
1059 res[1] = f.GetParameter(1);
1060
1061 return res;
1062 */
1063}
1064
1065// --------------------------------------------------------------------------
1066//
1067// Scale some image parameter plots using the scale factor and plot them
1068// together with the corresponding MC histograms.
1069// Called from DisplaySize
1070//
1071Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
1072{
1073 TString same(name);
1074 same += "Same";
1075
1076 TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
1077 TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
1078 if (!h1 || !h2)
1079 return kFALSE;
1080
1081 TObject *obj = plist.FindObject(plot);
1082 if (!obj)
1083 {
1084 *fLog << warn << plot << " not in parameter list... skipping." << endl;
1085 return kFALSE;
1086 }
1087
1088 TH1 *h3 = (TH1*)obj->FindObject(name);
1089 if (!h3)
1090 {
1091 *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
1092 return kFALSE;
1093 }
1094
1095
1096 const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
1097 const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
1098
1099 gPad->SetBorderMode(0);
1100 h2->SetLineColor(kBlack);
1101 h3->SetLineColor(kBlue);
1102 h2->Add(h1, -ascale);
1103
1104 //h2->Scale(1./ontime); //h2->Integral());
1105 h3->Scale(scale); //h3->Integral());
1106
1107 h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
1108
1109 h2 = h2->DrawCopy();
1110 h3 = h3->DrawCopy("same");
1111
1112 // Don't do this on the original object!
1113 h2->SetStats(kFALSE);
1114 h3->SetStats(kFALSE);
1115
1116 return kTRUE;
1117}
1118
1119// --------------------------------------------------------------------------
1120//
1121// Take a lot of histograms and plot them together in one plot.
1122// Calls PlotSame
1123//
1124Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
1125{
1126 *fLog << inf << "Reading from file: " << fPathIn << endl;
1127
1128 TFile file(fPathIn, "READ");
1129 if (!file.IsOpen())
1130 {
1131 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
1132 return kFALSE;
1133 }
1134
1135 file.cd();
1136 MStatusArray arr;
1137 if (arr.Read()<=0)
1138 {
1139 *fLog << "MStatusDisplay not found in file... abort." << endl;
1140 return kFALSE;
1141 }
1142
1143 TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
1144 if (!excess)
1145 return kFALSE;
1146
1147 // ------------------- Plot excess versus size -------------------
1148
1149 TCanvas &c = fDisplay->AddTab("Excess");
1150 c.Divide(3,2);
1151 c.cd(1);
1152 gPad->SetBorderMode(0);
1153 gPad->SetLogx();
1154 gPad->SetLogy();
1155 gPad->SetGridx();
1156 gPad->SetGridy();
1157
1158 excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
1159 excess = excess->DrawCopy("E2");
1160 // Don't do this on the original object!
1161 excess->SetStats(kFALSE);
1162 excess->SetMarkerStyle(kFullDotMedium);
1163 excess->SetFillColor(kBlack);
1164 excess->SetFillStyle(0);
1165 excess->SetName("Excess ");
1166 excess->SetDirectory(0);
1167
1168 TObject *o=0;
1169 if ((o=plist.FindObject("ExcessMC")))
1170 {
1171 TH1 *histsel = (TH1F*)o->FindObject("");
1172 if (histsel)
1173 {
1174 if (scale<0)
1175 scale = excess->Integral()/histsel->Integral();
1176
1177 histsel->Scale(scale);
1178 histsel->SetLineColor(kBlue);
1179 histsel->SetBit(kCanDelete);
1180 histsel = histsel->DrawCopy("E1 same");
1181 // Don't do this on the original object!
1182 histsel->SetStats(kFALSE);
1183
1184 fLog->Separator("Kolmogorov Test");
1185 histsel->KolmogorovTest(excess, "DX");
1186 fLog->Separator("Chi^2 Test");
1187 const Double_t p = histsel->Chi2Test(excess, "P");
1188
1189 TLatex tex;
1190 tex.SetBit(TLatex::kTextNDC);
1191 tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
1192 }
1193 }
1194
1195 // -------------- Comparison of Image Parameters --------------
1196 c.cd(2);
1197 PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
1198
1199 c.cd(3);
1200 PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
1201
1202 c.cd(4);
1203 PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
1204
1205 c.cd(5);
1206 PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
1207
1208 c.cd(6);
1209 PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
1210
1211 return kTRUE;
1212}
1213
1214// --------------------------------------------------------------------------
1215//
1216void MJSpectrum::DisplayCutEfficiency(const MHCollectionArea &area0, const MHCollectionArea &area1) const
1217{
1218 if (!fDisplay)
1219 return;
1220
1221 const TH1D &trig = area0.GetHEnergy();
1222 TH1D &cut = (TH1D&)*area1.GetHEnergy().Clone();
1223
1224 fDisplay->AddTab("CutEff");
1225
1226 gPad->SetBorderMode(0);
1227 gPad->SetFrameBorderMode(0);
1228 gPad->SetLogx();
1229 gPad->SetGridx();
1230 gPad->SetGridy();
1231
1232 cut.Divide(&trig);
1233 cut.Scale(100);
1234 cut.SetNameTitle("CutEff", "Background Supression: Cut efficiency (after star)");
1235 cut.SetYTitle("\\eta [%]");
1236 cut.SetDirectory(0);
1237 cut.SetMinimum(0);
1238 cut.SetMaximum(100);
1239 cut.SetBit(kCanDelete);
1240 cut.Draw();
1241
1242 TLine line;
1243 line.SetLineColor(kBlue);
1244 line.SetLineWidth(2);
1245 line.SetLineStyle(kDashed);
1246 line.DrawLine(cut.GetBinLowEdge(1), 50, cut.GetBinLowEdge(cut.GetNbinsX()+1), 50);
1247}
1248
1249void MJSpectrum::SetupHistEvtDist(MHn &hist) const
1250{
1251 hist.AddHist("MMcEvt.fEnergy");
1252 hist.InitName("EnergyDist;EnergyEst");
1253 hist.InitTitle("Unweighted event distribution (Real statistics);E [GeV];Counts;");
1254
1255 hist.AddHist("MPointingPos.fZd");
1256 hist.InitName("ZdDist;Theta");
1257 hist.InitTitle("Unweighted event distribution (Real statistics);Zd [\\circ];Counts;");
1258}
1259
1260void MJSpectrum::SetupHistEnergyEst(MHn &hist) const
1261{
1262 const char *res = "log10(MEnergyEst.fVal)-log10(MMcEvt.fEnergy)";
1263
1264 hist.AddHist("MHillas.fSize", res);
1265 hist.InitName("ResSize;Size;EnergyResidual");
1266 hist.InitTitle(";S [phe];\\Delta lg E;");
1267 hist.SetDrawOption("colz profx");
1268
1269 hist.AddHist("MPointingPos.fZd", res);
1270 hist.InitName("ResTheta;Theta;EnergyResidual");
1271 hist.InitTitle(";Zd [\\circ];\\Delta lg E;");
1272 hist.SetDrawOption("colz profx");
1273
1274 hist.AddHist("MNewImagePar.fLeakage1", res);
1275 hist.InitName("ResLeak;Leakage;EnergyResidual");
1276 hist.InitTitle(";Leak;\\Delta lg E;");
1277 hist.SetDrawOption("colz profx");
1278
1279 hist.AddHist("MHillasSrc.fDist*3.37e-3", res);
1280 hist.InitName("ResDist;Dist;EnergyResidual");
1281 hist.InitTitle(";D [\\circ];\\Delta lg E;");
1282 hist.SetDrawOption("colz profx");
1283}
1284
1285void MJSpectrum::SetupHistDisp(MHn &hist) const
1286{
1287 const char *res = "-Disp.fVal*sign(MHillasSrc.fCosDeltaAlpha)-MHillasSrc.fDist*3.37e-3";
1288
1289 hist.AddHist("MHillas.fSize", res);
1290 hist.InitName("ResSize;Size;ResidualDist");
1291 hist.InitTitle(";S [phe];Disp-Dist [\\circ];");
1292 hist.SetDrawOption("colz profx");
1293
1294 hist.AddHist("MPointingPos.fZd", res);
1295 hist.InitName("ResTheta;Theta;ResidualDist");
1296 hist.InitTitle(";Zd [\\circ];Disp-Dist [\\circ];");
1297 hist.SetDrawOption("colz profx");
1298
1299 hist.AddHist("MNewImagePar.fLeakage1", res);
1300 hist.InitName("ResLeak;Leakage;ResidualDist");
1301 hist.InitTitle(";Leak;Disp-Dist [\\circ];");
1302 hist.SetDrawOption("colz profx");
1303
1304 hist.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3", res);
1305 hist.InitName("ResSlope;Slope;ResidualDist");
1306 hist.InitTitle(";Slope;Disp-Dist [\\circ];");
1307 hist.SetDrawOption("colz profx");
1308}
1309
1310Bool_t MJSpectrum::Process(const MDataSet &set)
1311{
1312 if (!set.IsValid())
1313 {
1314 *fLog << err << "ERROR - DataSet invalid!" << endl;
1315 return kFALSE;
1316 }
1317
1318 if (!HasWritePermission(GetPathOut()))
1319 return kFALSE;
1320
1321 CheckEnv();
1322
1323 // --------------------------------------------------------------------------------
1324
1325 *fLog << inf;
1326 fLog->Separator(GetDescriptor());
1327 *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
1328 *fLog << endl;
1329
1330 if (fDisplay)
1331 fDisplay->SetWindowName(fName);
1332
1333 // Setup everything which is read from the ganymed file
1334 MBinning bins1("BinningAlpha");
1335 MBinning bins2("BinningEnergyEst");
1336 MBinning bins3("BinningTheta");
1337 MBinning bins4("BinningFalseSource");
1338 MBinning bins5("BinningWidth");
1339 MBinning bins6("BinningLength");
1340 MBinning bins7("BinningDist");
1341 MBinning bins8("BinningM3Long");
1342 MBinning bins9("BinningM3Trans");
1343 MBinning bins0("BinningSlope");
1344 MBinning binsa("BinningAsym");
1345 MBinning binsb("BinningConc1");
1346
1347 MAlphaFitter fit;
1348
1349 MParList plist;
1350 plist.AddToList(&bins0);
1351 plist.AddToList(&bins1);
1352 plist.AddToList(&bins3);
1353 plist.AddToList(&bins4);
1354 plist.AddToList(&bins5);
1355 plist.AddToList(&bins6);
1356 plist.AddToList(&bins7);
1357 plist.AddToList(&bins8);
1358 plist.AddToList(&bins9);
1359 plist.AddToList(&binsa);
1360 plist.AddToList(&binsb);
1361 plist.AddToList(&fit);
1362
1363 // Read from the ganymed file
1364 TH1D htheta, size;
1365 const Float_t ontime = ReadInput(plist, htheta, size);
1366 if (ontime<0)
1367 {
1368 *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
1369 return kFALSE;
1370 }
1371
1372 // Set Zenith angle binning to binning from the ganymed-file
1373 bins3.SetEdges(htheta, 'x');
1374
1375 // Read energy binning from resource file
1376 if (!CheckEnv(bins2))
1377 {
1378 *fLog << err << "ERROR - Reading energy binning from resources failed." << endl;
1379 return kFALSE;
1380 }
1381 plist.AddToList(&bins2); // For later use in MC processing
1382
1383 // Initialize weighting to a new spectrum as defined in the resource file
1384 MMcSpectrumWeight weight;
1385 if (!InitWeighting(set, weight))
1386 return kFALSE;
1387
1388 // Print Theta and energy binning for cross-checks
1389 *fLog << all << endl;
1390 bins2.Print();
1391 bins3.Print();
1392
1393 // Now we read the MC distribution as produced by corsika
1394 // vs zenith angle and energy.
1395 // Weight for the new energy spectrum defined in MMcSpectumWeight
1396 // are applied.
1397 // Also correction for different lower energy bounds and
1398 // different production areas (impact parameters) are applied.
1399 TH2D hist;
1400 hist.UseCurrentStyle();
1401 MH::SetBinning(&hist, &bins3, &bins2);
1402 if (!ReadOrigMCDistribution(set, hist, weight))
1403 return kFALSE;
1404
1405 // Check if user has closed the display
1406 if (!fDisplay)
1407 return kTRUE;
1408
1409 // Display histograms and calculate za-weights into htheta
1410 GetThetaDistribution(htheta, hist);
1411
1412 // Give the zenoith angle weights to the weighting task
1413 weight.SetWeightsZd(&htheta);
1414
1415 // No we apply the the zenith-angle-weights to the corsika produced
1416 // MC distribution. Unfortunately this must be done manually
1417 // because we are multiplying column by column
1418 for (int y=0; y<=hist.GetNbinsY()+1; y++)
1419 for (int x=0; x<=hist.GetNbinsX()+1; x++)
1420 {
1421 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*htheta.GetBinContent(x));
1422 hist.SetBinError(x, y, hist.GetBinError(x, y) *htheta.GetBinContent(x));
1423 }
1424
1425 // Display the resulting distribution and check it matches
1426 // the observation time distribution (this could be false
1427 // for example if you miss MCs of some zenith angles, which you have
1428 // data for)
1429 if (!DisplayResult(hist))
1430 return kFALSE;
1431
1432 // -------------- Fill excess events versus energy ---------------
1433
1434 // Refill excess histogram to determin the excess events
1435 TH1D excess;
1436 if (!Refill(plist, excess))
1437 return kFALSE;
1438
1439 // Print the setup and result of the MAlphaFitter, print used cuts
1440 PrintSetup(fit);
1441
1442 // ------------------------- Final loop --------------------------
1443
1444 *fLog << endl;
1445 fLog->Separator("Calculate efficiencies");
1446 *fLog << endl;
1447
1448 MTaskList tlist2;
1449 plist.Replace(&tlist2);
1450
1451 MReadMarsFile read("Events");
1452 read.DisableAutoScheme();
1453 if (!set.AddFilesOn(read))
1454 return kFALSE;
1455
1456 // Selector to get correct (final) theta-distribution
1457 //temp1.SetXTitle("MPointingPos.fZd");
1458 //
1459 //MH3 mh3(temp1);
1460 //
1461 //MFEventSelector2 sel2(mh3);
1462 //sel2.SetHistIsProbability();
1463 //
1464 //MContinue contsel(&sel2);
1465 //contsel.SetInverted();
1466
1467 // Get correct source position
1468 //MSrcPosCalc calc;
1469
1470 // Calculate corresponding Hillas parameters
1471 //MHillasCalc hcalc1;
1472 //MHillasCalc hcalc2("MHillasCalcAnti");
1473 //hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1474 //hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1475 //hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1476 //hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1477
1478 // Fill collection area and energy estimator (unfolding)
1479 // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1480 MHCollectionArea area0("TriggerArea");
1481 MHCollectionArea area1;
1482 area0.SetHistAll(hist);
1483 area1.SetHistAll(hist);
1484
1485 MHEnergyEst hest;
1486
1487 MFillH fill30(&area0, "", "FillTriggerArea");
1488 MFillH fill31(&area1, "", "FillCollectionArea");
1489 MFillH fill4(&hest, "", "FillEnergyEst");
1490 MFillH fill5("MHThreshold", "", "FillThreshold");
1491 fill30.SetWeight();
1492 fill31.SetWeight();
1493 fill4.SetWeight();
1494 fill5.SetWeight();
1495 fill30.SetNameTab("TrigArea");
1496 fill31.SetNameTab("ColArea");
1497 fill4.SetNameTab("E-Est");
1498 fill5.SetNameTab("Threshold");
1499
1500/*
1501 MH3 henergy("MMcEvt.fEnergy");
1502 henergy.SetName("EventDist;EnergyEst");
1503 henergy.SetTitle("Unweighted event distribution (Real statistics);E [GeV];Counts;");
1504 henergy.Sumw2();
1505 */
1506
1507 // ---------------------------------------------------------
1508
1509 MBinning binsA(50, 10, 100000, "BinningSize", "log");
1510 MBinning binsC(50, 0, 0.3, "BinningLeakage", "lin");
1511 MBinning binsB(51, -1, 1, "BinningEnergyResidual", "lin");
1512 MBinning binsD(51, -1, 1, "BinningResidualDist", "lin");
1513
1514 plist.AddToList(&binsA);
1515 plist.AddToList(&binsB);
1516 plist.AddToList(&binsC);
1517 plist.AddToList(&binsD);
1518
1519 MHn heest("Energy", "Energy Residual (lg E_{est} - lg E_{mc})");
1520 SetupHistEnergyEst(heest);
1521
1522 //MHn hdisp("Disp", "Dist residual (Disp-Dist)");
1523 //SetupHistDisp(hdisp);
1524
1525 MHn henergy("EvtDist");
1526 SetupHistEvtDist(henergy);
1527
1528 MFillH fill4b(&heest, "", "FillEnergyResidual");
1529 fill4b.SetWeight();
1530
1531 //MFillH fill4c(&hdisp, "", "FillDispResidual");
1532 //fill4c.SetWeight();
1533
1534 // ---------------------------------------------------------
1535
1536 MH3 hsize("MHillas.fSize");
1537 hsize.SetName("ExcessMC");
1538 hsize.Sumw2();
1539
1540 MBinning bins(size, "BinningExcessMC");
1541 plist.AddToList(&hsize);
1542 plist.AddToList(&bins);
1543
1544 // ---------------------------------------------------------
1545
1546 MFillH fill0a(&henergy, "", "FillEventDist");
1547 MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1548 MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1549 MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1550 MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1551 MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1552 MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1553 MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1554 MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1555 fill0a.SetNameTab("EvtDist");
1556 fill1a.SetNameTab("PreCut");
1557 fill2a.SetNameTab("PostCut");
1558 fill3a.SetNameTab("VsSize");
1559 fill4a.SetNameTab("HilExt");
1560 fill5a.SetNameTab("HilSrc");
1561 fill6a.SetNameTab("ImgPar");
1562 fill7a.SetNameTab("NewPar");
1563 fill8a.SetBit(MFillH::kDoNotDisplay);
1564 fill1a.SetWeight();
1565 fill2a.SetWeight();
1566 fill3a.SetWeight();
1567 fill4a.SetWeight();
1568 fill5a.SetWeight();
1569 fill6a.SetWeight();
1570 fill7a.SetWeight();
1571 fill8a.SetWeight();
1572
1573 // FIXME: To be done: A task checking the lower 1% after the lower
1574 // energy limit!
1575
1576 MTaskEnv taskenv0("CalcDisp");
1577 taskenv0.SetDefault(fCalcDisp);
1578
1579 MTaskEnv taskenv1("CalcHadronness");
1580 taskenv1.SetDefault(fCalcHadronness);
1581
1582 MEnergyEstimate est;
1583 MTaskEnv taskenv2("EstimateEnergy");
1584 taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1585
1586 tlist2.AddToList(&read);
1587 // If no weighting should be applied but the events should
1588 // be thrown away according to the theta distribution
1589 // it is enabled here
1590 //if (!fRawMc && fNoThetaWeights)
1591 // tlist2.AddToList(&contsel);
1592 //tlist2.AddToList(&calc);
1593 //tlist2.AddToList(&hcalc1);
1594 //tlist2.AddToList(&hcalc2);
1595 tlist2.AddToList(&weight);
1596 tlist2.AddToList(&fill1a);
1597 tlist2.AddToList(&fill30);
1598 tlist2.AddToList(fCutQ);
1599 tlist2.AddToList(fCut0);
1600 tlist2.AddToList(&taskenv0);
1601 tlist2.AddToList(&taskenv1);
1602 tlist2.AddToList(fCut1);
1603 tlist2.AddToList(fCutS);
1604 tlist2.AddToList(fCut2);
1605 tlist2.AddToList(fCut3);
1606 tlist2.AddToList(&taskenv2);
1607 tlist2.AddToList(&fill31);
1608 tlist2.AddToList(&fill4);
1609 tlist2.AddToList(&fill4b);
1610 tlist2.AddToList(&fill5);
1611 tlist2.AddToList(&fill0a);
1612 tlist2.AddToList(&fill2a);
1613 tlist2.AddToList(&fill3a);
1614 tlist2.AddToList(&fill4a);
1615 tlist2.AddToList(&fill5a);
1616 tlist2.AddToList(&fill6a);
1617 tlist2.AddToList(&fill7a);
1618 tlist2.AddToList(&fill8a);
1619 //tlist2.AddToList(&fill9a);
1620
1621 // by setting it here it is distributed to all consecutive tasks
1622 tlist2.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
1623
1624 MEvtLoop loop2(fName); // ***** fName *****
1625 loop2.SetParList(&plist);
1626 loop2.SetDisplay(fDisplay);
1627 loop2.SetLogStream(fLog);
1628
1629 if (!SetupEnv(loop2))
1630 return kFALSE;
1631
1632 if (!loop2.Eventloop(fMaxEvents))
1633 {
1634 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1635 return kFALSE;
1636 }
1637
1638 if (!loop2.GetDisplay())
1639 {
1640 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1641 return kFALSE;
1642 }
1643
1644 gLog.Separator("Energy Estimator");
1645 if (plist.FindObject("EstimateEnergy"))
1646 plist.FindObject("EstimateEnergy")->Print();
1647
1648 gLog.Separator("Spectrum");
1649
1650 // -------------------------- Spectrum ----------------------------
1651
1652 // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1653 TArrayD res(DisplaySpectrum(area1, excess, hest, ontime));
1654
1655 // Spectrum fitted (convert res[1] from TeV to GeV)
1656 TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
1657
1658 // Number of events this spectrum would produce per s and m^2
1659 Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1660
1661 // scale with effective collection area to get the event rate (N/s)
1662 // scale with the effective observation time to absolute observed events
1663 n *= area1.GetCollectionAreaAbs()*ontime; // N
1664
1665 // Now calculate the scale factor from the number of events
1666 // produced and the number of events which should have been
1667 // observed with our telescope in the time ontime
1668 const Double_t scale = n/area1.GetEntries();
1669
1670 // Print normalization constant
1671 cout << "MC normalization factor: " << scale << endl;
1672
1673 // Display cut efficiency
1674 DisplayCutEfficiency(area0, area1);
1675
1676 // Overlay normalized plots
1677 DisplaySize(plist, scale);
1678
1679 // check if output should be written
1680 if (fPathOut.IsNull())
1681 return kTRUE;
1682
1683 // Write the output
1684 TObjArray cont;
1685 cont.Add((TObject*)GetEnv());
1686 cont.Add(&area0);
1687 cont.Add(&area1);
1688 cont.Add(&hest);
1689
1690 if (fDisplay)
1691 cont.Add(fDisplay);
1692
1693 if (!WriteContainer(cont, "", "RECREATE"))
1694 {
1695 *fLog << err << GetDescriptor() << ": Writing result failed." << endl;
1696 return kFALSE;
1697 }
1698
1699 *fLog << all << GetDescriptor() << ": Done." << endl;
1700 *fLog << endl << endl;
1701
1702 return kTRUE;
1703}
Note: See TracBrowser for help on using the repository browser.