1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appear in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | !
|
---|
18 | ! Author(s): Thomas Bretz, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2005
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | /////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MJSpectrum
|
---|
28 | //
|
---|
29 | // Program to calculate spectrum
|
---|
30 | //
|
---|
31 | /////////////////////////////////////////////////////////////////////////////
|
---|
32 | #include "MJSpectrum.h"
|
---|
33 |
|
---|
34 | // Root
|
---|
35 | #include <TF1.h>
|
---|
36 | #include <TH1.h>
|
---|
37 | #include <TH2.h>
|
---|
38 | #include <TFile.h>
|
---|
39 | #include <TChain.h>
|
---|
40 | #include <TCanvas.h>
|
---|
41 | #include <TObjArray.h>
|
---|
42 |
|
---|
43 | // Environment
|
---|
44 | #include "MLog.h"
|
---|
45 | #include "MLogManip.h"
|
---|
46 |
|
---|
47 | #include "MStatusArray.h"
|
---|
48 | #include "MStatusDisplay.h"
|
---|
49 |
|
---|
50 | // Container
|
---|
51 | #include "MH3.h"
|
---|
52 | #include "MBinning.h"
|
---|
53 | #include "MDataSet.h"
|
---|
54 |
|
---|
55 | // Spectrum
|
---|
56 | #include "../mhflux/MAlphaFitter.h"
|
---|
57 | #include "../mhflux/MHAlpha.h"
|
---|
58 | #include "../mhflux/MHCollectionArea.h"
|
---|
59 | #include "../mhflux/MHEnergyEst.h"
|
---|
60 |
|
---|
61 | // Eventloop
|
---|
62 | #include "MEvtLoop.h"
|
---|
63 | #include "MTaskList.h"
|
---|
64 | #include "MParList.h"
|
---|
65 |
|
---|
66 | // Tasks/Filter
|
---|
67 | #include "MReadMarsFile.h"
|
---|
68 | #include "MReadMarsFile.h"
|
---|
69 | #include "MFEventSelector2.h"
|
---|
70 | #include "MFDataMember.h"
|
---|
71 | #include "MEnergyEstimate.h"
|
---|
72 | #include "MTaskEnv.h"
|
---|
73 | #include "MFillH.h"
|
---|
74 | #include "MHillasCalc.h"
|
---|
75 | #include "MSrcPosCalc.h"
|
---|
76 | #include "MContinue.h"
|
---|
77 |
|
---|
78 | ClassImp(MJSpectrum);
|
---|
79 |
|
---|
80 | using namespace std;
|
---|
81 |
|
---|
82 | MJSpectrum::MJSpectrum(const char *name, const char *title)
|
---|
83 | : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
|
---|
84 | fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
|
---|
85 | {
|
---|
86 | fName = name ? name : "MJSpectrum";
|
---|
87 | fTitle = title ? title : "Standard program to calculate spectrum";
|
---|
88 | }
|
---|
89 |
|
---|
90 | MJSpectrum::~MJSpectrum()
|
---|
91 | {
|
---|
92 | if (fCut0)
|
---|
93 | delete fCut0;
|
---|
94 | if (fCut1)
|
---|
95 | delete fCut1;
|
---|
96 | if (fCut2)
|
---|
97 | delete fCut2;
|
---|
98 | if (fCut3)
|
---|
99 | delete fCut3;
|
---|
100 | if (fEstimateEnergy)
|
---|
101 | delete fEstimateEnergy;
|
---|
102 | }
|
---|
103 |
|
---|
104 | // --------------------------------------------------------------------------
|
---|
105 | //
|
---|
106 | // Setup a task estimating the energy. The given task is cloned.
|
---|
107 | //
|
---|
108 | void MJSpectrum::SetEnergyEstimator(const MTask *task)
|
---|
109 | {
|
---|
110 | if (fEstimateEnergy)
|
---|
111 | delete fEstimateEnergy;
|
---|
112 | fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
|
---|
113 | }
|
---|
114 |
|
---|
115 | Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
|
---|
116 | {
|
---|
117 | if (task)
|
---|
118 | {
|
---|
119 | delete task;
|
---|
120 | task = 0;
|
---|
121 | }
|
---|
122 |
|
---|
123 | task = (MTask*)gFile->Get(name);
|
---|
124 | if (!task)
|
---|
125 | {
|
---|
126 | *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
|
---|
127 | return kFALSE;
|
---|
128 | }
|
---|
129 | if (!task->InheritsFrom(MTask::Class()))
|
---|
130 | {
|
---|
131 | *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
|
---|
132 | delete task;
|
---|
133 | return kFALSE;
|
---|
134 | }
|
---|
135 |
|
---|
136 | task->SetName(name);
|
---|
137 |
|
---|
138 | return kTRUE;
|
---|
139 | }
|
---|
140 |
|
---|
141 | void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
|
---|
142 | {
|
---|
143 | fLog->Separator("Alpha Fitter");
|
---|
144 | *fLog << all;
|
---|
145 | fit.Print();
|
---|
146 |
|
---|
147 | fLog->Separator("Used Cuts");
|
---|
148 | fCut0->Print();
|
---|
149 | fCut1->Print();
|
---|
150 | fCut2->Print();
|
---|
151 | fCut3->Print();
|
---|
152 |
|
---|
153 | //gLog.Separator("Energy Estimator");
|
---|
154 | //fEstimateEnergy->Print();
|
---|
155 | }
|
---|
156 |
|
---|
157 | Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
|
---|
158 | {
|
---|
159 | *fLog << inf << "Reading from file: " << fPathIn << endl;
|
---|
160 |
|
---|
161 | TFile file(fPathIn, "READ");
|
---|
162 | if (!file.IsOpen())
|
---|
163 | {
|
---|
164 | *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
|
---|
165 | return -1;
|
---|
166 | }
|
---|
167 |
|
---|
168 | MStatusArray arr;
|
---|
169 | if (arr.Read()<=0)
|
---|
170 | {
|
---|
171 | *fLog << "MStatusDisplay not found in file... abort." << endl;
|
---|
172 | return -1;
|
---|
173 | }
|
---|
174 |
|
---|
175 | TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
|
---|
176 | TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
|
---|
177 | if (!vstime || !size)
|
---|
178 | return -1;
|
---|
179 |
|
---|
180 | vstime->Copy(h1);
|
---|
181 | size->Copy(h2);
|
---|
182 | h1.SetDirectory(0);
|
---|
183 | h2.SetDirectory(0);
|
---|
184 |
|
---|
185 | if (fDisplay)
|
---|
186 | arr.DisplayIn(*fDisplay, "Hist");
|
---|
187 |
|
---|
188 | if (!ReadTask(fCut0, "Cut0"))
|
---|
189 | return -1;
|
---|
190 | if (!ReadTask(fCut1, "Cut1"))
|
---|
191 | return -1;
|
---|
192 | if (!ReadTask(fCut2, "Cut2"))
|
---|
193 | return -1;
|
---|
194 | if (!ReadTask(fCut3, "Cut3"))
|
---|
195 | return -1;
|
---|
196 |
|
---|
197 | TObjArray arrread;
|
---|
198 |
|
---|
199 | TIter Next(plist);
|
---|
200 | TObject *o=0;
|
---|
201 | while ((o=Next()))
|
---|
202 | if (o->InheritsFrom(MBinning::Class()))
|
---|
203 | arrread.Add(o);
|
---|
204 |
|
---|
205 | arrread.Add(plist.FindObject("MAlphaFitter"));
|
---|
206 |
|
---|
207 | if (!ReadContainer(arrread))
|
---|
208 | return -1;
|
---|
209 |
|
---|
210 | return vstime->Integral();
|
---|
211 | }
|
---|
212 |
|
---|
213 | Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
|
---|
214 | {
|
---|
215 | // Some debug output
|
---|
216 | fLog->Separator("Compiling original MC distribution");
|
---|
217 |
|
---|
218 | *fLog << inf << "Please stand by, this may take a while..." << flush;
|
---|
219 |
|
---|
220 | if (fDisplay)
|
---|
221 | fDisplay->SetStatusLine1("Compiling MC distribution...");
|
---|
222 |
|
---|
223 | // Create chain
|
---|
224 | TChain chain("OriginalMC");
|
---|
225 | set.AddFilesOn(chain);
|
---|
226 |
|
---|
227 | // Prepare histogram
|
---|
228 | h.Reset();
|
---|
229 |
|
---|
230 | // Fill histogram from chain
|
---|
231 | h.SetDirectory(gROOT);
|
---|
232 | if (h.InheritsFrom(TH2::Class()))
|
---|
233 | {
|
---|
234 | h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
|
---|
235 | h.SetXTitle("\\Theta [\\circ]");
|
---|
236 | h.SetYTitle("E [GeV]");
|
---|
237 | h.SetZTitle("Counts");
|
---|
238 | chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
|
---|
239 | }
|
---|
240 | else
|
---|
241 | {
|
---|
242 | h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
|
---|
243 | h.SetXTitle("\\Theta [\\circ]");
|
---|
244 | h.SetYTitle("Counts");
|
---|
245 | chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
|
---|
246 | }
|
---|
247 | h.SetDirectory(0);
|
---|
248 |
|
---|
249 | *fLog << "done." << endl;
|
---|
250 | if (fDisplay)
|
---|
251 | fDisplay->SetStatusLine2("done.");
|
---|
252 |
|
---|
253 | if (h.GetEntries()>0)
|
---|
254 | return kTRUE;
|
---|
255 |
|
---|
256 | *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
|
---|
257 |
|
---|
258 | return h.GetEntries()>0;
|
---|
259 | }
|
---|
260 |
|
---|
261 | Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
|
---|
262 | {
|
---|
263 | // Display some stuff
|
---|
264 | if (fDisplay)
|
---|
265 | {
|
---|
266 | TCanvas &c = fDisplay->AddTab("ZdDist");
|
---|
267 | c.Divide(2,2);
|
---|
268 |
|
---|
269 | // On-Time vs. Theta
|
---|
270 | c.cd(1);
|
---|
271 | gPad->SetBorderMode(0);
|
---|
272 | temp1.DrawCopy();
|
---|
273 |
|
---|
274 | // Number of MC events (produced) vs Theta
|
---|
275 | c.cd(2);
|
---|
276 | gPad->SetBorderMode(0);
|
---|
277 | temp2.SetName("NVsTheta");
|
---|
278 | temp2.DrawCopy();
|
---|
279 |
|
---|
280 | c.cd(4);
|
---|
281 | gPad->SetBorderMode(0);
|
---|
282 |
|
---|
283 | c.cd(3);
|
---|
284 | gPad->SetBorderMode(0);
|
---|
285 | }
|
---|
286 |
|
---|
287 | // Calculate the Probability
|
---|
288 | temp1.Divide(&temp2);
|
---|
289 | temp1.Scale(1./temp1.GetMaximum());
|
---|
290 |
|
---|
291 | // Some cosmetics: Name, Axis, etc.
|
---|
292 | temp1.SetName("ProbVsTheta");
|
---|
293 | temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
|
---|
294 | temp1.SetYTitle("Probability");
|
---|
295 | if (fDisplay)
|
---|
296 | temp1.DrawCopy();
|
---|
297 |
|
---|
298 | return kTRUE;
|
---|
299 | }
|
---|
300 |
|
---|
301 | void MJSpectrum::DisplayResult(const TH2D &h2) const
|
---|
302 | {
|
---|
303 | if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
|
---|
304 | return;
|
---|
305 |
|
---|
306 | TH1D &proj = *h2.ProjectionX();
|
---|
307 | proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
|
---|
308 | proj.SetXTitle("\\Theta [\\circ]");
|
---|
309 | proj.SetYTitle("Counts");
|
---|
310 | proj.SetLineColor(kBlue);
|
---|
311 | proj.SetDirectory(0);
|
---|
312 | proj.SetBit(kCanDelete);
|
---|
313 |
|
---|
314 | TVirtualPad *pad = gPad;
|
---|
315 |
|
---|
316 | pad->cd(4);
|
---|
317 | proj.DrawCopy();
|
---|
318 |
|
---|
319 | pad->cd(1);
|
---|
320 | TH1D *theta = (TH1D*)gPad->FindObject("Theta");
|
---|
321 | if (theta)
|
---|
322 | {
|
---|
323 | proj.Scale(theta->GetMaximum()/proj.GetMaximum());
|
---|
324 | theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
|
---|
325 | }
|
---|
326 | proj.Draw("same");
|
---|
327 | }
|
---|
328 |
|
---|
329 | Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
|
---|
330 | {
|
---|
331 | // Try to find the class used to determin the signal!
|
---|
332 | TString cls("MHAlpha");
|
---|
333 | if (fDisplay)
|
---|
334 | {
|
---|
335 | TCanvas *c = fDisplay->GetCanvas("Hist");
|
---|
336 | if (c)
|
---|
337 | {
|
---|
338 | TIter Next(c->GetListOfPrimitives());
|
---|
339 | TObject *obj=0;
|
---|
340 | while ((obj=Next()))
|
---|
341 | if (obj->InheritsFrom(MHAlpha::Class()))
|
---|
342 | break;
|
---|
343 | if (obj)
|
---|
344 | cls = obj->ClassName();
|
---|
345 | }
|
---|
346 | }
|
---|
347 |
|
---|
348 | cout << "FOUND: "<< cls << endl;
|
---|
349 |
|
---|
350 | // Now fill the histogram
|
---|
351 | *fLog << endl;
|
---|
352 | fLog->Separator("Refill Excess");
|
---|
353 | *fLog << endl;
|
---|
354 |
|
---|
355 | MTaskList tlist;
|
---|
356 | plist.AddToList(&tlist);
|
---|
357 |
|
---|
358 | MReadTree read("Events");
|
---|
359 | read.DisableAutoScheme();
|
---|
360 | read.AddFile(fPathIn);
|
---|
361 |
|
---|
362 | MEnergyEstimate est;
|
---|
363 | MTaskEnv taskenv1("EstimateEnergy");
|
---|
364 | taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
|
---|
365 |
|
---|
366 | MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
|
---|
367 | MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
|
---|
368 |
|
---|
369 | MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
|
---|
370 | MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
|
---|
371 |
|
---|
372 | fill1.SetFilter(&f0);
|
---|
373 | fill2.SetFilter(&f1);
|
---|
374 |
|
---|
375 | tlist.AddToList(&read);
|
---|
376 | tlist.AddToList(&taskenv1);
|
---|
377 | tlist.AddToList(&f0);
|
---|
378 | tlist.AddToList(&f1);
|
---|
379 | tlist.AddToList(&fill1);
|
---|
380 | tlist.AddToList(&fill2);
|
---|
381 |
|
---|
382 | MEvtLoop loop(fName);
|
---|
383 | loop.SetParList(&plist);
|
---|
384 | loop.SetDisplay(fDisplay);
|
---|
385 | loop.SetLogStream(fLog);
|
---|
386 |
|
---|
387 | if (!SetupEnv(loop))
|
---|
388 | return kFALSE;
|
---|
389 |
|
---|
390 | if (!loop.Eventloop())
|
---|
391 | {
|
---|
392 | *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
|
---|
393 | return kFALSE;
|
---|
394 | }
|
---|
395 |
|
---|
396 | if (!loop.GetDisplay())
|
---|
397 | {
|
---|
398 | *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
|
---|
399 | return kFALSE;
|
---|
400 | }
|
---|
401 |
|
---|
402 | const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
|
---|
403 | if (!halpha)
|
---|
404 | {
|
---|
405 | *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
|
---|
406 | return kFALSE;
|
---|
407 | }
|
---|
408 |
|
---|
409 | halpha->GetHEnergy().Copy(h2);
|
---|
410 | h2.SetDirectory(0);
|
---|
411 |
|
---|
412 | return kTRUE;
|
---|
413 | }
|
---|
414 |
|
---|
415 | Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
|
---|
416 | {
|
---|
417 | MTaskList tlist1;
|
---|
418 | plist.Replace(&tlist1);
|
---|
419 |
|
---|
420 | MReadMarsFile readmc("OriginalMC");
|
---|
421 | //readmc.DisableAutoScheme();
|
---|
422 | set.AddFilesOn(readmc);
|
---|
423 | readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
|
---|
424 | readmc.EnableBranch("MMcEvtBasic.fEnergy");
|
---|
425 |
|
---|
426 | mh1.SetLogy();
|
---|
427 | mh1.SetLogz();
|
---|
428 | mh1.SetName("ThetaE");
|
---|
429 |
|
---|
430 | MFillH fill0(&mh1);
|
---|
431 | //fill0.SetDrawOption("projx only");
|
---|
432 |
|
---|
433 | MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
|
---|
434 | MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
|
---|
435 | if (bins2 && bins3)
|
---|
436 | {
|
---|
437 | bins2->SetName("BinningThetaEY");
|
---|
438 | bins3->SetName("BinningThetaEX");
|
---|
439 | }
|
---|
440 | tlist1.AddToList(&readmc);
|
---|
441 |
|
---|
442 | temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
|
---|
443 | MH3 mh3mc(temp1);
|
---|
444 |
|
---|
445 | MFEventSelector2 sel1(mh3mc);
|
---|
446 | sel1.SetHistIsProbability();
|
---|
447 |
|
---|
448 | fill0.SetFilter(&sel1);
|
---|
449 |
|
---|
450 | if (!fRawMc)
|
---|
451 | tlist1.AddToList(&sel1);
|
---|
452 | tlist1.AddToList(&fill0);
|
---|
453 |
|
---|
454 | MEvtLoop loop1(fName);
|
---|
455 | loop1.SetParList(&plist);
|
---|
456 | loop1.SetLogStream(fLog);
|
---|
457 | loop1.SetDisplay(fDisplay);
|
---|
458 |
|
---|
459 | if (!SetupEnv(loop1))
|
---|
460 | return kFALSE;
|
---|
461 |
|
---|
462 | if (!loop1.Eventloop(fMaxEvents))
|
---|
463 | {
|
---|
464 | *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
|
---|
465 | return kFALSE;
|
---|
466 | }
|
---|
467 |
|
---|
468 | if (!loop1.GetDisplay())
|
---|
469 | {
|
---|
470 | *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
|
---|
471 | return kFALSE;
|
---|
472 | }
|
---|
473 |
|
---|
474 | if (bins2 && bins3)
|
---|
475 | {
|
---|
476 | bins2->SetName("BinningEnergyEst");
|
---|
477 | bins3->SetName("BinningTheta");
|
---|
478 | }
|
---|
479 | return kTRUE;
|
---|
480 | }
|
---|
481 |
|
---|
482 | void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
|
---|
483 | {
|
---|
484 | TH1D collarea(area.GetHEnergy());
|
---|
485 | TH1D spectrum(excess);
|
---|
486 | TH1D weights;
|
---|
487 | hest.GetWeights(weights);
|
---|
488 |
|
---|
489 | cout << "Effective On time: " << ontime << "s" << endl;
|
---|
490 |
|
---|
491 | spectrum.SetDirectory(NULL);
|
---|
492 | spectrum.SetBit(kCanDelete);
|
---|
493 | spectrum.Scale(1./ontime);
|
---|
494 | spectrum.Divide(&collarea);
|
---|
495 | spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
|
---|
496 | spectrum.SetYTitle("N/sm^{2}");
|
---|
497 |
|
---|
498 | TCanvas &c1 = fDisplay->AddTab("Spectrum");
|
---|
499 | c1.Divide(2,2);
|
---|
500 | c1.cd(1);
|
---|
501 | gPad->SetBorderMode(0);
|
---|
502 | gPad->SetLogx();
|
---|
503 | gPad->SetLogy();
|
---|
504 | gPad->SetGridx();
|
---|
505 | gPad->SetGridy();
|
---|
506 | collarea.DrawCopy();
|
---|
507 |
|
---|
508 | c1.cd(2);
|
---|
509 | gPad->SetBorderMode(0);
|
---|
510 | gPad->SetLogx();
|
---|
511 | gPad->SetLogy();
|
---|
512 | gPad->SetGridx();
|
---|
513 | gPad->SetGridy();
|
---|
514 | spectrum.DrawCopy();
|
---|
515 |
|
---|
516 | c1.cd(3);
|
---|
517 | gPad->SetBorderMode(0);
|
---|
518 | gPad->SetLogx();
|
---|
519 | gPad->SetLogy();
|
---|
520 | gPad->SetGridx();
|
---|
521 | gPad->SetGridy();
|
---|
522 | weights.DrawCopy();
|
---|
523 |
|
---|
524 | //spectrum.Divide(&weights);
|
---|
525 | spectrum.Multiply(&weights);
|
---|
526 | spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
|
---|
527 |
|
---|
528 | for (int i=0; i<excess.GetNbinsX(); i++)
|
---|
529 | {
|
---|
530 | spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
|
---|
531 | spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
|
---|
532 | }
|
---|
533 |
|
---|
534 | c1.cd(4);
|
---|
535 | gPad->SetBorderMode(0);
|
---|
536 | gPad->SetLogx();
|
---|
537 | gPad->SetLogy();
|
---|
538 | gPad->SetGridx();
|
---|
539 | gPad->SetGridy();
|
---|
540 | spectrum.SetXTitle("E [GeV]");
|
---|
541 | spectrum.SetYTitle("N/TeVsm^{2}");
|
---|
542 | spectrum.DrawCopy();
|
---|
543 |
|
---|
544 | TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
|
---|
545 | f.SetParameter(0, -2.87);
|
---|
546 | f.SetParameter(1, 1.9e-6);
|
---|
547 | f.SetLineColor(kGreen);
|
---|
548 | spectrum.Fit(&f, "NI", "", 55, 2e4);
|
---|
549 | f.DrawCopy("same");
|
---|
550 |
|
---|
551 | /*
|
---|
552 | TString str;
|
---|
553 | str += "(1.68#pm0.15)10^{-7}";
|
---|
554 | str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
|
---|
555 | str += "\\frac{ph}{TeVm^{2}s}";
|
---|
556 |
|
---|
557 | TLatex tex;
|
---|
558 | tex.DrawLatex(2e2, 7e-5, str);
|
---|
559 | */
|
---|
560 | }
|
---|
561 |
|
---|
562 | Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
|
---|
563 | {
|
---|
564 | cout << name << endl;
|
---|
565 | TString same(name);
|
---|
566 | same += "Same";
|
---|
567 |
|
---|
568 | TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
|
---|
569 | TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
|
---|
570 | if (!h1 || !h2)
|
---|
571 | return kFALSE;
|
---|
572 |
|
---|
573 | TObject *obj = plist.FindObject(plot);
|
---|
574 | if (!obj)
|
---|
575 | {
|
---|
576 | *fLog << warn << plot << " not in parameter list... skipping." << endl;
|
---|
577 | return kFALSE;
|
---|
578 | }
|
---|
579 |
|
---|
580 | TH1 *h3 = (TH1*)obj->FindObject(name);
|
---|
581 | if (!h3)
|
---|
582 | {
|
---|
583 | *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
|
---|
584 | return kFALSE;
|
---|
585 | }
|
---|
586 |
|
---|
587 |
|
---|
588 | const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
|
---|
589 | const Double_t scale = fit ? fit->GetScaleFactor() : 1;
|
---|
590 |
|
---|
591 | gPad->SetBorderMode(0);
|
---|
592 | h2->SetLineColor(kBlack);
|
---|
593 | h3->SetLineColor(kBlue);
|
---|
594 | h2->Add(h1, -scale);
|
---|
595 |
|
---|
596 | h2->Scale(1./h2->Integral());
|
---|
597 | h3->Scale(1./h3->Integral());
|
---|
598 |
|
---|
599 | h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
|
---|
600 |
|
---|
601 | h2 = h2->DrawCopy();
|
---|
602 | h3 = h3->DrawCopy("same");
|
---|
603 |
|
---|
604 | // Don't do this on the original object!
|
---|
605 | h2->SetStats(kFALSE);
|
---|
606 | h3->SetStats(kFALSE);
|
---|
607 |
|
---|
608 | return kTRUE;
|
---|
609 | }
|
---|
610 |
|
---|
611 | Bool_t MJSpectrum::DisplaySize(MParList &plist) const
|
---|
612 | {
|
---|
613 | *fLog << inf << "Reading from file: " << fPathIn << endl;
|
---|
614 |
|
---|
615 | TFile file(fPathIn, "READ");
|
---|
616 | if (!file.IsOpen())
|
---|
617 | {
|
---|
618 | *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
|
---|
619 | return kFALSE;
|
---|
620 | }
|
---|
621 |
|
---|
622 | file.cd();
|
---|
623 | MStatusArray arr;
|
---|
624 | if (arr.Read()<=0)
|
---|
625 | {
|
---|
626 | *fLog << "MStatusDisplay not found in file... abort." << endl;
|
---|
627 | return kFALSE;
|
---|
628 | }
|
---|
629 |
|
---|
630 | TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
|
---|
631 | if (!excess)
|
---|
632 | return kFALSE;
|
---|
633 |
|
---|
634 | // ------------------- Plot excess versus size -------------------
|
---|
635 |
|
---|
636 | TCanvas &c = fDisplay->AddTab("Excess");
|
---|
637 | c.Divide(3,2);
|
---|
638 | c.cd(1);
|
---|
639 | gPad->SetBorderMode(0);
|
---|
640 | gPad->SetLogx();
|
---|
641 | gPad->SetLogy();
|
---|
642 | gPad->SetGridx();
|
---|
643 | gPad->SetGridy();
|
---|
644 |
|
---|
645 | excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
|
---|
646 | excess->Scale(1./excess->Integral());
|
---|
647 | excess = excess->DrawCopy();
|
---|
648 | // Don't do this on the original object!
|
---|
649 | excess->SetStats(kFALSE);
|
---|
650 |
|
---|
651 | TObject *o=0;
|
---|
652 | if ((o=plist.FindObject("ExcessSize")))
|
---|
653 | {
|
---|
654 | TH1 *histsel = (TH1F*)o->FindObject("");
|
---|
655 | if (histsel)
|
---|
656 | {
|
---|
657 | histsel->Scale(1./histsel->Integral());
|
---|
658 | histsel->SetLineColor(kBlue);
|
---|
659 | histsel->SetBit(kCanDelete);
|
---|
660 | histsel = histsel->DrawCopy("same");
|
---|
661 | // Don't do this on the original object!
|
---|
662 | histsel->SetStats(kFALSE);
|
---|
663 | }
|
---|
664 | }
|
---|
665 |
|
---|
666 | // -------------- Comparison of Image Parameters --------------
|
---|
667 | c.cd(2);
|
---|
668 | PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost");
|
---|
669 |
|
---|
670 | c.cd(3);
|
---|
671 | PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
|
---|
672 |
|
---|
673 | c.cd(4);
|
---|
674 | PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost");
|
---|
675 |
|
---|
676 | c.cd(5);
|
---|
677 | PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost");
|
---|
678 |
|
---|
679 | c.cd(6);
|
---|
680 | PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost");
|
---|
681 |
|
---|
682 | return kTRUE;
|
---|
683 | }
|
---|
684 |
|
---|
685 | Bool_t MJSpectrum::Process(const MDataSet &set)
|
---|
686 | {
|
---|
687 | if (!set.IsValid())
|
---|
688 | {
|
---|
689 | *fLog << err << "ERROR - DataSet invalid!" << endl;
|
---|
690 | return kFALSE;
|
---|
691 | }
|
---|
692 |
|
---|
693 | CheckEnv();
|
---|
694 |
|
---|
695 | // --------------------------------------------------------------------------------
|
---|
696 |
|
---|
697 | *fLog << inf;
|
---|
698 | fLog->Separator(GetDescriptor());
|
---|
699 | *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
|
---|
700 | *fLog << endl;
|
---|
701 |
|
---|
702 | MBinning bins1("BinningAlpha");
|
---|
703 | MBinning bins2("BinningEnergyEst");
|
---|
704 | MBinning bins3("BinningTheta");
|
---|
705 | MBinning bins4("BinningFalseSource");
|
---|
706 | MBinning bins5("BinningWidth");
|
---|
707 | MBinning bins6("BinningLength");
|
---|
708 | MBinning bins7("BinningDist");
|
---|
709 | MBinning bins8("BinningMaxDist");
|
---|
710 | MBinning bins9("BinningM3Long");
|
---|
711 | MBinning bins0("BinningConc1");
|
---|
712 |
|
---|
713 | MAlphaFitter fit;
|
---|
714 |
|
---|
715 | MParList plist;
|
---|
716 | plist.AddToList(&bins1);
|
---|
717 | plist.AddToList(&bins2);
|
---|
718 | plist.AddToList(&bins3);
|
---|
719 | plist.AddToList(&bins4);
|
---|
720 | plist.AddToList(&bins5);
|
---|
721 | plist.AddToList(&bins6);
|
---|
722 | plist.AddToList(&bins7);
|
---|
723 | plist.AddToList(&bins8);
|
---|
724 | plist.AddToList(&bins9);
|
---|
725 | plist.AddToList(&bins0);
|
---|
726 | plist.AddToList(&fit);
|
---|
727 |
|
---|
728 | TH1D temp1, size;
|
---|
729 | const Float_t ontime = ReadInput(plist, temp1, size);
|
---|
730 | if (ontime<0)
|
---|
731 | {
|
---|
732 | *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
|
---|
733 | return kFALSE;
|
---|
734 | }
|
---|
735 |
|
---|
736 | PrintSetup(fit);
|
---|
737 | bins3.SetEdges(temp1, 'x');
|
---|
738 |
|
---|
739 | TH1D temp2(temp1);
|
---|
740 | if (!ReadOrigMCDistribution(set, temp2))
|
---|
741 | return kFALSE;
|
---|
742 |
|
---|
743 | if (!GetThetaDistribution(temp1, temp2))
|
---|
744 | return kFALSE;
|
---|
745 |
|
---|
746 | TH1D excess;
|
---|
747 | if (!Refill(plist, excess))
|
---|
748 | return kFALSE;
|
---|
749 |
|
---|
750 | TH2D hist;
|
---|
751 | MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
|
---|
752 | if (fSimpleMode)
|
---|
753 | {
|
---|
754 | hist.UseCurrentStyle();
|
---|
755 | MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
|
---|
756 | if (!ReadOrigMCDistribution(set, hist))
|
---|
757 | return kFALSE;
|
---|
758 |
|
---|
759 | if (!fRawMc)
|
---|
760 | {
|
---|
761 | for (int y=0; y<hist.GetNbinsY(); y++)
|
---|
762 | for (int x=0; x<hist.GetNbinsX(); x++)
|
---|
763 | hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
|
---|
764 | }
|
---|
765 | }
|
---|
766 | else
|
---|
767 | if (!IntermediateLoop(plist, mh1, temp1, set))
|
---|
768 | return kFALSE;
|
---|
769 |
|
---|
770 | DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
|
---|
771 |
|
---|
772 | // ------------------------- Final loop --------------------------
|
---|
773 |
|
---|
774 | *fLog << endl;
|
---|
775 | fLog->Separator("Calculate efficiencies");
|
---|
776 | *fLog << endl;
|
---|
777 |
|
---|
778 | MTaskList tlist2;
|
---|
779 | plist.Replace(&tlist2);
|
---|
780 |
|
---|
781 | MReadMarsFile read("Events");
|
---|
782 | read.DisableAutoScheme();
|
---|
783 | set.AddFilesOn(read);
|
---|
784 |
|
---|
785 | // Selector to get correct (final) theta-distribution
|
---|
786 | temp1.SetXTitle("MPointingPos.fZd");
|
---|
787 |
|
---|
788 | MH3 mh3(temp1);
|
---|
789 |
|
---|
790 | MFEventSelector2 sel2(mh3);
|
---|
791 | sel2.SetHistIsProbability();
|
---|
792 |
|
---|
793 | MContinue contsel(&sel2);
|
---|
794 | contsel.SetInverted();
|
---|
795 |
|
---|
796 | // Get correct source position
|
---|
797 | MSrcPosCalc calc;
|
---|
798 |
|
---|
799 | // Calculate corresponding Hillas parameters
|
---|
800 | MHillasCalc hcalc1;
|
---|
801 | MHillasCalc hcalc2("MHillasCalcAnti");
|
---|
802 | hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
|
---|
803 | hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
|
---|
804 | hcalc2.SetNameHillasSrc("MHillasSrcAnti");
|
---|
805 | hcalc2.SetNameSrcPosCam("MSrcPosAnti");
|
---|
806 |
|
---|
807 | // Fill collection area and energy estimator (unfolding)
|
---|
808 | // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
|
---|
809 | MHCollectionArea area;
|
---|
810 | area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
|
---|
811 | MHEnergyEst hest;
|
---|
812 |
|
---|
813 | MFillH fill3(&area, "", "FillCollectionArea");
|
---|
814 | MFillH fill4(&hest, "", "FillEnergyEst");
|
---|
815 |
|
---|
816 | MH3 hsize("MHillas.fSize");
|
---|
817 | //MH3 henergy("MEnergyEst.fVal");
|
---|
818 | hsize.SetName("ExcessSize");
|
---|
819 | //henergy.SetName("EnergyEst");
|
---|
820 | MBinning bins(size, "BinningExcessSize");
|
---|
821 | plist.AddToList(&hsize);
|
---|
822 | //plist.AddToList(&henergy);
|
---|
823 | plist.AddToList(&bins);
|
---|
824 |
|
---|
825 | MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
|
---|
826 | MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
|
---|
827 | MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
|
---|
828 | MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
|
---|
829 | MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
|
---|
830 | MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
|
---|
831 | MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
|
---|
832 | MFillH fill8a("ExcessSize [MH3]", "", "FillExcessSize");
|
---|
833 | //MFillH fill9a("EnergyEst [MH3]", "", "FillExcessEEst");
|
---|
834 | fill1a.SetNameTab("PreCut");
|
---|
835 | fill2a.SetNameTab("PostCut");
|
---|
836 | fill3a.SetNameTab("VsSize");
|
---|
837 | fill4a.SetNameTab("HilExt");
|
---|
838 | fill5a.SetNameTab("HilSrc");
|
---|
839 | fill6a.SetNameTab("ImgPar");
|
---|
840 | fill7a.SetNameTab("NewPar");
|
---|
841 | fill8a.SetBit(MFillH::kDoNotDisplay);
|
---|
842 | //fill9a.SetBit(MFillH::kDoNotDisplay);
|
---|
843 |
|
---|
844 | MEnergyEstimate est;
|
---|
845 | MTaskEnv taskenv1("EstimateEnergy");
|
---|
846 | taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
|
---|
847 |
|
---|
848 | tlist2.AddToList(&read);
|
---|
849 | if (!fRawMc)
|
---|
850 | tlist2.AddToList(&contsel);
|
---|
851 | tlist2.AddToList(&calc);
|
---|
852 | tlist2.AddToList(&hcalc1);
|
---|
853 | tlist2.AddToList(&hcalc2);
|
---|
854 | tlist2.AddToList(&fill1a);
|
---|
855 | tlist2.AddToList(fCut0);
|
---|
856 | tlist2.AddToList(fCut1);
|
---|
857 | tlist2.AddToList(fCut2);
|
---|
858 | tlist2.AddToList(fCut3);
|
---|
859 | tlist2.AddToList(&taskenv1);
|
---|
860 | tlist2.AddToList(&fill3);
|
---|
861 | tlist2.AddToList(&fill4);
|
---|
862 | tlist2.AddToList(&fill2a);
|
---|
863 | tlist2.AddToList(&fill3a);
|
---|
864 | tlist2.AddToList(&fill4a);
|
---|
865 | tlist2.AddToList(&fill5a);
|
---|
866 | tlist2.AddToList(&fill6a);
|
---|
867 | tlist2.AddToList(&fill7a);
|
---|
868 | tlist2.AddToList(&fill8a);
|
---|
869 | //tlist2.AddToList(&fill9a);
|
---|
870 |
|
---|
871 | MEvtLoop loop2(fName);
|
---|
872 | loop2.SetParList(&plist);
|
---|
873 | loop2.SetDisplay(fDisplay);
|
---|
874 | loop2.SetLogStream(fLog);
|
---|
875 |
|
---|
876 | if (!SetupEnv(loop2))
|
---|
877 | return kFALSE;
|
---|
878 |
|
---|
879 | if (!loop2.Eventloop(fMaxEvents))
|
---|
880 | {
|
---|
881 | *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
|
---|
882 | return kFALSE;
|
---|
883 | }
|
---|
884 |
|
---|
885 | if (!loop2.GetDisplay())
|
---|
886 | {
|
---|
887 | *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
|
---|
888 | return kFALSE;
|
---|
889 | }
|
---|
890 |
|
---|
891 | // -------------------------- Spectrum ----------------------------
|
---|
892 |
|
---|
893 | DisplaySpectrum(area, excess, hest, ontime);
|
---|
894 | DisplaySize(plist);
|
---|
895 |
|
---|
896 | if (!fPathOut.IsNull())
|
---|
897 | fDisplay->SaveAsRoot(fPathOut);
|
---|
898 |
|
---|
899 | return kTRUE;
|
---|
900 | }
|
---|