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

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