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

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