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

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