| 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 "MFEventSelector2.h"
|
|---|
| 69 | #include "MFDataMember.h"
|
|---|
| 70 | #include "MFillH.h"
|
|---|
| 71 | #include "MHillasCalc.h"
|
|---|
| 72 | #include "MSrcPosCalc.h"
|
|---|
| 73 | #include "MContinue.h"
|
|---|
| 74 |
|
|---|
| 75 | ClassImp(MJSpectrum);
|
|---|
| 76 |
|
|---|
| 77 | using namespace std;
|
|---|
| 78 |
|
|---|
| 79 | MJSpectrum::MJSpectrum(const char *name, const char *title)
|
|---|
| 80 | : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
|
|---|
| 81 | fRefill(kFALSE), fSimpleMode(kTRUE)
|
|---|
| 82 | {
|
|---|
| 83 | fName = name ? name : "MJSpectrum";
|
|---|
| 84 | fTitle = title ? title : "Standard program to calculate spectrum";
|
|---|
| 85 | }
|
|---|
| 86 |
|
|---|
| 87 | MJSpectrum::~MJSpectrum()
|
|---|
| 88 | {
|
|---|
| 89 | if (fCut0)
|
|---|
| 90 | delete fCut0;
|
|---|
| 91 | if (fCut1)
|
|---|
| 92 | delete fCut1;
|
|---|
| 93 | if (fCut2)
|
|---|
| 94 | delete fCut2;
|
|---|
| 95 | if (fCut3)
|
|---|
| 96 | delete fCut3;
|
|---|
| 97 | if (fEstimateEnergy)
|
|---|
| 98 | delete fEstimateEnergy;
|
|---|
| 99 | }
|
|---|
| 100 |
|
|---|
| 101 | Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
|
|---|
| 102 | {
|
|---|
| 103 | if (task)
|
|---|
| 104 | {
|
|---|
| 105 | delete task;
|
|---|
| 106 | task = 0;
|
|---|
| 107 | }
|
|---|
| 108 |
|
|---|
| 109 | task = (MTask*)gFile->Get(name);
|
|---|
| 110 | if (!task)
|
|---|
| 111 | {
|
|---|
| 112 | *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
|
|---|
| 113 | return kFALSE;
|
|---|
| 114 | }
|
|---|
| 115 | if (!task->InheritsFrom(MTask::Class()))
|
|---|
| 116 | {
|
|---|
| 117 | *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
|
|---|
| 118 | delete task;
|
|---|
| 119 | return kFALSE;
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | task->SetName(name);
|
|---|
| 123 |
|
|---|
| 124 | return kTRUE;
|
|---|
| 125 | }
|
|---|
| 126 |
|
|---|
| 127 | Bool_t MJSpectrum::ReadInput(const MParList &plist)
|
|---|
| 128 | {
|
|---|
| 129 | const TString fname = fPathIn;
|
|---|
| 130 |
|
|---|
| 131 | *fLog << inf << "Reading from file: " << fname << endl;
|
|---|
| 132 |
|
|---|
| 133 | TFile file(fname, "READ");
|
|---|
| 134 | if (!file.IsOpen())
|
|---|
| 135 | {
|
|---|
| 136 | *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
|
|---|
| 137 | return kFALSE;
|
|---|
| 138 | }
|
|---|
| 139 |
|
|---|
| 140 | if (!ReadTask(fCut0, "Cut0"))
|
|---|
| 141 | return kFALSE;
|
|---|
| 142 | if (!ReadTask(fCut1, "Cut1"))
|
|---|
| 143 | return kFALSE;
|
|---|
| 144 | if (!ReadTask(fCut2, "Cut2"))
|
|---|
| 145 | return kFALSE;
|
|---|
| 146 | if (!ReadTask(fCut3, "Cut3"))
|
|---|
| 147 | return kFALSE;
|
|---|
| 148 |
|
|---|
| 149 | if (!ReadTask(fEstimateEnergy, "EstimateEnergy"))
|
|---|
| 150 | return kFALSE;
|
|---|
| 151 |
|
|---|
| 152 | TObjArray arr;
|
|---|
| 153 |
|
|---|
| 154 | TIter Next(plist);
|
|---|
| 155 | TObject *o=0;
|
|---|
| 156 | while ((o=Next()))
|
|---|
| 157 | if (o->InheritsFrom(MBinning::Class()))
|
|---|
| 158 | arr.Add(o);
|
|---|
| 159 |
|
|---|
| 160 | arr.Add(plist.FindObject("MAlphaFitter"));
|
|---|
| 161 |
|
|---|
| 162 | return ReadContainer(arr);
|
|---|
| 163 | }
|
|---|
| 164 |
|
|---|
| 165 | void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
|
|---|
| 166 | {
|
|---|
| 167 | gLog.Separator("Alpha Fitter");
|
|---|
| 168 | *fLog << all;
|
|---|
| 169 | fit.Print();
|
|---|
| 170 |
|
|---|
| 171 | gLog.Separator("Used Cuts");
|
|---|
| 172 | fCut0->Print();
|
|---|
| 173 | fCut1->Print();
|
|---|
| 174 | fCut2->Print();
|
|---|
| 175 | fCut3->Print();
|
|---|
| 176 |
|
|---|
| 177 | gLog.Separator("Energy Estimator");
|
|---|
| 178 | fEstimateEnergy->Print();
|
|---|
| 179 | }
|
|---|
| 180 |
|
|---|
| 181 | Float_t MJSpectrum::ReadHistograms(TH1D &h1, TH1D &h2) const
|
|---|
| 182 | {
|
|---|
| 183 | TFile file(fPathIn, "READ");
|
|---|
| 184 | if (!file.IsOpen())
|
|---|
| 185 | {
|
|---|
| 186 | *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
|
|---|
| 187 | return -1;
|
|---|
| 188 | }
|
|---|
| 189 |
|
|---|
| 190 | MStatusArray *arr = (MStatusArray*)file.Get("MStatusDisplay");
|
|---|
| 191 | if (!arr)
|
|---|
| 192 | {
|
|---|
| 193 | gLog << "MStatusDisplay not found in file... abort." << endl;
|
|---|
| 194 | return -1;
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | TH1D *vstime = (TH1D*)arr->FindObjectInCanvas("Theta", "TH1D", "OnTime");
|
|---|
| 198 | TH1D *excen = (TH1D*)arr->FindObjectInCanvas("ExcessEnergy", "TH1D", "MHAlpha");
|
|---|
| 199 | if (!vstime || !excen)
|
|---|
| 200 | return -1;
|
|---|
| 201 |
|
|---|
| 202 | vstime->Copy(h1);
|
|---|
| 203 | excen->Copy(h2);
|
|---|
| 204 | h1.SetDirectory(0);
|
|---|
| 205 | h2.SetDirectory(0);
|
|---|
| 206 |
|
|---|
| 207 | if (fDisplay)
|
|---|
| 208 | arr->DisplayIn(*fDisplay, "MHAlpha");
|
|---|
| 209 |
|
|---|
| 210 | delete arr;
|
|---|
| 211 |
|
|---|
| 212 | return vstime->Integral();
|
|---|
| 213 | }
|
|---|
| 214 |
|
|---|
| 215 | Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
|
|---|
| 216 | {
|
|---|
| 217 | // Some debug output
|
|---|
| 218 | fLog->Separator("Compiling original MC distribution");
|
|---|
| 219 |
|
|---|
| 220 | *fLog << inf << "Please stand by, this may take a while..." << flush;
|
|---|
| 221 |
|
|---|
| 222 | if (fDisplay)
|
|---|
| 223 | fDisplay->SetStatusLine1("Compiling MC distribution...");
|
|---|
| 224 |
|
|---|
| 225 | // Create chain
|
|---|
| 226 | TChain chain("OriginalMC");
|
|---|
| 227 | set.AddFilesOn(chain);
|
|---|
| 228 |
|
|---|
| 229 | // Prepare histogram
|
|---|
| 230 | h.Reset();
|
|---|
| 231 |
|
|---|
| 232 | // Fill histogram from chain
|
|---|
| 233 | h.SetDirectory(gROOT);
|
|---|
| 234 | if (h.InheritsFrom(TH2::Class()))
|
|---|
| 235 | {
|
|---|
| 236 | h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
|
|---|
| 237 | h.SetXTitle("\\Theta [\\circ]");
|
|---|
| 238 | h.SetYTitle("E [GeV]");
|
|---|
| 239 | h.SetZTitle("Counts");
|
|---|
| 240 | chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
|
|---|
| 241 | }
|
|---|
| 242 | else
|
|---|
| 243 | {
|
|---|
| 244 | h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
|
|---|
| 245 | h.SetXTitle("\\Theta [\\circ]");
|
|---|
| 246 | h.SetYTitle("Counts");
|
|---|
| 247 | chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
|
|---|
| 248 | }
|
|---|
| 249 | h.SetDirectory(0);
|
|---|
| 250 |
|
|---|
| 251 | *fLog << "done." << endl;
|
|---|
| 252 | if (fDisplay)
|
|---|
| 253 | fDisplay->SetStatusLine2("done.");
|
|---|
| 254 |
|
|---|
| 255 | if (h.GetEntries()>0)
|
|---|
| 256 | return kTRUE;
|
|---|
| 257 |
|
|---|
| 258 | *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
|
|---|
| 259 |
|
|---|
| 260 | return h.GetEntries()>0;
|
|---|
| 261 | }
|
|---|
| 262 |
|
|---|
| 263 | Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
|
|---|
| 264 | {
|
|---|
| 265 | // Display some stuff
|
|---|
| 266 | if (fDisplay)
|
|---|
| 267 | {
|
|---|
| 268 | TCanvas &c = fDisplay->AddTab("ZdDist");
|
|---|
| 269 | c.Divide(2,2);
|
|---|
| 270 |
|
|---|
| 271 | // On-Time vs. Theta
|
|---|
| 272 | c.cd(1);
|
|---|
| 273 | gPad->SetBorderMode(0);
|
|---|
| 274 | temp1.DrawCopy();
|
|---|
| 275 |
|
|---|
| 276 | // Number of MC events (produced) vs Theta
|
|---|
| 277 | c.cd(2);
|
|---|
| 278 | gPad->SetBorderMode(0);
|
|---|
| 279 | temp2.SetName("NVsTheta");
|
|---|
| 280 | temp2.DrawCopy();
|
|---|
| 281 |
|
|---|
| 282 | c.cd(4);
|
|---|
| 283 | gPad->SetBorderMode(0);
|
|---|
| 284 |
|
|---|
| 285 | c.cd(3);
|
|---|
| 286 | gPad->SetBorderMode(0);
|
|---|
| 287 | }
|
|---|
| 288 |
|
|---|
| 289 | // Calculate the Probability
|
|---|
| 290 | temp1.Divide(&temp2);
|
|---|
| 291 | temp1.Scale(1./temp1.GetMaximum());
|
|---|
| 292 |
|
|---|
| 293 | // Some cosmetics: Name, Axis, etc.
|
|---|
| 294 | temp1.SetName("ProbVsTheta");
|
|---|
| 295 | temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
|
|---|
| 296 | temp1.SetYTitle("Probability");
|
|---|
| 297 | if (fDisplay)
|
|---|
| 298 | temp1.DrawCopy();
|
|---|
| 299 |
|
|---|
| 300 | return kTRUE;
|
|---|
| 301 | }
|
|---|
| 302 |
|
|---|
| 303 | void MJSpectrum::DisplayResult(const MH3 &mh1) const
|
|---|
| 304 | {
|
|---|
| 305 | if (!fDisplay->CdCanvas("ZdDist"))
|
|---|
| 306 | return;
|
|---|
| 307 |
|
|---|
| 308 | const TH2D &h2 = (const TH2D&)mh1.GetHist();
|
|---|
| 309 |
|
|---|
| 310 | TH1D &proj = *h2.ProjectionX();
|
|---|
| 311 | proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
|
|---|
| 312 | proj.SetXTitle("\\Theta [\\circ]");
|
|---|
| 313 | proj.SetYTitle("Counts");
|
|---|
| 314 | proj.SetLineColor(kBlue);
|
|---|
| 315 | proj.SetDirectory(0);
|
|---|
| 316 | proj.SetBit(kCanDelete);
|
|---|
| 317 |
|
|---|
| 318 | TVirtualPad *pad = gPad;
|
|---|
| 319 |
|
|---|
| 320 | pad->cd(4);
|
|---|
| 321 | proj.DrawCopy();
|
|---|
| 322 |
|
|---|
| 323 | pad->cd(1);
|
|---|
| 324 | TH1D *theta = (TH1D*)gPad->FindObject("Theta");
|
|---|
| 325 | if (theta)
|
|---|
| 326 | {
|
|---|
| 327 | proj.Scale(theta->Integral()/proj.Integral());
|
|---|
| 328 | theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
|
|---|
| 329 | }
|
|---|
| 330 | proj.Draw("same");
|
|---|
| 331 | }
|
|---|
| 332 |
|
|---|
| 333 | Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
|
|---|
| 334 | {
|
|---|
| 335 | *fLog << endl;
|
|---|
| 336 | fLog->Separator("Refill Excess");
|
|---|
| 337 | *fLog << endl;
|
|---|
| 338 |
|
|---|
| 339 | MTaskList tlist;
|
|---|
| 340 | plist.AddToList(&tlist);
|
|---|
| 341 |
|
|---|
| 342 | MReadTree read("Events");
|
|---|
| 343 | read.DisableAutoScheme();
|
|---|
| 344 | read.AddFile(fPathIn);
|
|---|
| 345 |
|
|---|
| 346 | MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
|
|---|
| 347 | MFillH fill2("MHAlpha", "MHillasSrc", "FillAlpha");
|
|---|
| 348 |
|
|---|
| 349 | MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
|
|---|
| 350 | MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
|
|---|
| 351 |
|
|---|
| 352 | fill1.SetFilter(&f0);
|
|---|
| 353 | fill2.SetFilter(&f1);
|
|---|
| 354 |
|
|---|
| 355 | tlist.AddToList(&read);
|
|---|
| 356 | tlist.AddToList(fEstimateEnergy);
|
|---|
| 357 | tlist.AddToList(&f0);
|
|---|
| 358 | tlist.AddToList(&f1);
|
|---|
| 359 | tlist.AddToList(&fill1);
|
|---|
| 360 | tlist.AddToList(&fill2);
|
|---|
| 361 |
|
|---|
| 362 | MEvtLoop loop(fName);
|
|---|
| 363 | loop.SetParList(&plist);
|
|---|
| 364 | loop.SetDisplay(fDisplay);
|
|---|
| 365 | loop.SetLogStream(fLog);
|
|---|
| 366 |
|
|---|
| 367 | if (!SetupEnv(loop))
|
|---|
| 368 | return kFALSE;
|
|---|
| 369 |
|
|---|
| 370 | if (!loop.Eventloop())
|
|---|
| 371 | {
|
|---|
| 372 | *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
|
|---|
| 373 | return kFALSE;
|
|---|
| 374 | }
|
|---|
| 375 |
|
|---|
| 376 | tlist.PrintStatistics();
|
|---|
| 377 |
|
|---|
| 378 | if (!loop.GetDisplay())
|
|---|
| 379 | {
|
|---|
| 380 | *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
|
|---|
| 381 | return kFALSE;
|
|---|
| 382 | }
|
|---|
| 383 |
|
|---|
| 384 | const MHAlpha *halpha = (MHAlpha *)plist.FindObject("MHAlpha");
|
|---|
| 385 | if (!halpha)
|
|---|
| 386 | {
|
|---|
| 387 | *fLog << err << GetDescriptor() << ": MHAlpha not found... abort." << endl;
|
|---|
| 388 | return kFALSE;
|
|---|
| 389 | }
|
|---|
| 390 |
|
|---|
| 391 | halpha->GetHEnergy().Copy(h2);
|
|---|
| 392 | h2.SetDirectory(0);
|
|---|
| 393 |
|
|---|
| 394 | return kTRUE;
|
|---|
| 395 | }
|
|---|
| 396 |
|
|---|
| 397 | Bool_t MJSpectrum::Process(const MDataSet &set)
|
|---|
| 398 | {
|
|---|
| 399 | if (!set.IsValid())
|
|---|
| 400 | {
|
|---|
| 401 | *fLog << err << "ERROR - DataSet invalid!" << endl;
|
|---|
| 402 | return kFALSE;
|
|---|
| 403 | }
|
|---|
| 404 |
|
|---|
| 405 | CheckEnv();
|
|---|
| 406 |
|
|---|
| 407 | // --------------------------------------------------------------------------------
|
|---|
| 408 |
|
|---|
| 409 | *fLog << inf;
|
|---|
| 410 | fLog->Separator(GetDescriptor());
|
|---|
| 411 | *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
|
|---|
| 412 | *fLog << endl;
|
|---|
| 413 |
|
|---|
| 414 | MBinning bins1("BinningAlpha");
|
|---|
| 415 | MBinning bins2("BinningEnergyEst");
|
|---|
| 416 | MBinning bins3("BinningTheta");
|
|---|
| 417 | MBinning bins4("BinningFalseSource");
|
|---|
| 418 | MBinning bins5("BinningWidth");
|
|---|
| 419 | MBinning bins6("BinningLength");
|
|---|
| 420 | MBinning bins7("BinningDist");
|
|---|
| 421 | MBinning bins8("BinningMaxDist");
|
|---|
| 422 |
|
|---|
| 423 | MAlphaFitter fit;
|
|---|
| 424 |
|
|---|
| 425 | MParList plist;
|
|---|
| 426 | plist.AddToList(&bins1);
|
|---|
| 427 | plist.AddToList(&bins2);
|
|---|
| 428 | plist.AddToList(&bins3);
|
|---|
| 429 | plist.AddToList(&bins4);
|
|---|
| 430 | plist.AddToList(&bins5);
|
|---|
| 431 | plist.AddToList(&bins6);
|
|---|
| 432 | plist.AddToList(&bins7);
|
|---|
| 433 | plist.AddToList(&bins8);
|
|---|
| 434 | plist.AddToList(&fit);
|
|---|
| 435 |
|
|---|
| 436 | if (!ReadInput(plist))
|
|---|
| 437 | return kFALSE;
|
|---|
| 438 |
|
|---|
| 439 | PrintSetup(fit);
|
|---|
| 440 |
|
|---|
| 441 | TH1D temp1, excess;
|
|---|
| 442 | const Float_t ontime = ReadHistograms(temp1, excess);
|
|---|
| 443 | if (ontime<0)
|
|---|
| 444 | {
|
|---|
| 445 | *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
|
|---|
| 446 | return kFALSE;
|
|---|
| 447 | }
|
|---|
| 448 | /*
|
|---|
| 449 | TH2D hist;
|
|---|
| 450 | if (fSimpleMode)
|
|---|
| 451 | {
|
|---|
| 452 | hist.UseCurrentStyle();
|
|---|
| 453 | MH::SetBinning(&hist, temp1.GetXaxis(), excess.GetXaxis());
|
|---|
| 454 | if (!ReadOrigMCDistribution(set, hist))
|
|---|
| 455 | return kFALSE;
|
|---|
| 456 |
|
|---|
| 457 | TCanvas &c =fDisplay->AddTab("test");
|
|---|
| 458 | c.Divide(2,2);
|
|---|
| 459 | c.cd(1);
|
|---|
| 460 | gPad->SetLogy();
|
|---|
| 461 | hist.DrawCopy();
|
|---|
| 462 | c.cd(2);
|
|---|
| 463 | hist.ProjectionX()->Draw();
|
|---|
| 464 | c.cd(3);
|
|---|
| 465 | hist.ProjectionY()->Draw();
|
|---|
| 466 |
|
|---|
| 467 | MH3 mh1(hist);
|
|---|
| 468 | mh1.SetLogy();
|
|---|
| 469 | mh1.SetLogz();
|
|---|
| 470 | mh1.SetName("ThetaE");
|
|---|
| 471 | }
|
|---|
| 472 | */
|
|---|
| 473 | TH1D temp2(temp1);
|
|---|
| 474 | if (!ReadOrigMCDistribution(set, temp2))
|
|---|
| 475 | return kFALSE;
|
|---|
| 476 |
|
|---|
| 477 | if (!GetThetaDistribution(temp1, temp2))
|
|---|
| 478 | return kFALSE;
|
|---|
| 479 |
|
|---|
| 480 | if (fRefill && !Refill(plist, temp2))
|
|---|
| 481 | return kFALSE;
|
|---|
| 482 |
|
|---|
| 483 | // ------------------- Loop1 ----------------------------
|
|---|
| 484 |
|
|---|
| 485 | MTaskList tlist1;
|
|---|
| 486 | plist.Replace(&tlist1);
|
|---|
| 487 |
|
|---|
| 488 | MReadMarsFile readmc("OriginalMC");
|
|---|
| 489 | //readmc.DisableAutoScheme();
|
|---|
| 490 | set.AddFilesOn(readmc);
|
|---|
| 491 | readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
|
|---|
| 492 | readmc.EnableBranch("MMcEvtBasic.fEnergy");
|
|---|
| 493 |
|
|---|
| 494 | MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
|
|---|
| 495 | mh1.SetLogy();
|
|---|
| 496 | mh1.SetLogz();
|
|---|
| 497 | mh1.SetName("ThetaE");
|
|---|
| 498 |
|
|---|
| 499 | MFillH fill0(&mh1);
|
|---|
| 500 | //fill0.SetDrawOption("projx only");
|
|---|
| 501 |
|
|---|
| 502 | MBinning binsx(bins3, "BinningThetaEX");
|
|---|
| 503 | MBinning binsy(bins2, "BinningThetaEY");
|
|---|
| 504 | plist.AddToList(&binsx);
|
|---|
| 505 | plist.AddToList(&binsy);
|
|---|
| 506 | tlist1.AddToList(&readmc);
|
|---|
| 507 |
|
|---|
| 508 | temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
|
|---|
| 509 | MH3 mh3mc(temp1);
|
|---|
| 510 |
|
|---|
| 511 | MFEventSelector2 sel1(mh3mc);
|
|---|
| 512 | sel1.SetHistIsProbability();
|
|---|
| 513 |
|
|---|
| 514 | fill0.SetFilter(&sel1);
|
|---|
| 515 |
|
|---|
| 516 | tlist1.AddToList(&sel1);
|
|---|
| 517 | tlist1.AddToList(&fill0);
|
|---|
| 518 |
|
|---|
| 519 | MEvtLoop loop1(fName);
|
|---|
| 520 | loop1.SetParList(&plist);
|
|---|
| 521 | loop1.SetLogStream(fLog);
|
|---|
| 522 | loop1.SetDisplay(fDisplay);
|
|---|
| 523 |
|
|---|
| 524 | if (!SetupEnv(loop1))
|
|---|
| 525 | return kFALSE;
|
|---|
| 526 |
|
|---|
| 527 | if (!loop1.Eventloop(fMaxEvents))
|
|---|
| 528 | {
|
|---|
| 529 | *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
|
|---|
| 530 | return kFALSE;
|
|---|
| 531 | }
|
|---|
| 532 |
|
|---|
| 533 | tlist1.PrintStatistics();
|
|---|
| 534 |
|
|---|
| 535 | if (!loop1.GetDisplay())
|
|---|
| 536 | {
|
|---|
| 537 | *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
|
|---|
| 538 | return kFALSE;
|
|---|
| 539 | }
|
|---|
| 540 |
|
|---|
| 541 | DisplayResult(mh1);
|
|---|
| 542 |
|
|---|
| 543 | // ------------------------- Final loop --------------------------
|
|---|
| 544 |
|
|---|
| 545 | *fLog << endl;
|
|---|
| 546 | fLog->Separator("Calculate efficiencies");
|
|---|
| 547 | *fLog << endl;
|
|---|
| 548 |
|
|---|
| 549 | MTaskList tlist2;
|
|---|
| 550 | plist.Replace(&tlist2);
|
|---|
| 551 |
|
|---|
| 552 | MReadMarsFile read("Events");
|
|---|
| 553 | read.DisableAutoScheme();
|
|---|
| 554 | set.AddFilesOn(read);
|
|---|
| 555 |
|
|---|
| 556 | // Selector to get correct (final) theta-distribution
|
|---|
| 557 | temp1.SetXTitle("MPointingPos.fZd");
|
|---|
| 558 |
|
|---|
| 559 | MH3 mh3(temp1);
|
|---|
| 560 |
|
|---|
| 561 | MFEventSelector2 sel2(mh3);
|
|---|
| 562 | sel2.SetHistIsProbability();
|
|---|
| 563 |
|
|---|
| 564 | MContinue contsel(&sel2);
|
|---|
| 565 | contsel.SetInverted();
|
|---|
| 566 |
|
|---|
| 567 | // Get correct source position
|
|---|
| 568 | MSrcPosCalc calc;
|
|---|
| 569 |
|
|---|
| 570 | // Calculate corresponding Hillas parameters
|
|---|
| 571 | MHillasCalc hcalc1;
|
|---|
| 572 | MHillasCalc hcalc2("MHillasCalcAnti");
|
|---|
| 573 | hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
|
|---|
| 574 | hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
|
|---|
| 575 | hcalc2.SetNameHillasSrc("MHillasSrcAnti");
|
|---|
| 576 | hcalc2.SetNameSrcPosCam("MSrcPosAnti");
|
|---|
| 577 |
|
|---|
| 578 | // Fill collection area and energy estimator (unfolding)
|
|---|
| 579 | // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
|
|---|
| 580 | MHCollectionArea area;
|
|---|
| 581 | area.SetHistAll((TH2D&)mh1.GetHist());
|
|---|
| 582 | MHEnergyEst hest;
|
|---|
| 583 |
|
|---|
| 584 | MFillH fill3(&area, "", "FillCollectionArea");
|
|---|
| 585 | MFillH fill4(&hest, "", "FillEnergyEst");
|
|---|
| 586 |
|
|---|
| 587 | MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
|
|---|
| 588 | MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
|
|---|
| 589 | MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
|
|---|
| 590 | MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
|
|---|
| 591 | MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
|
|---|
| 592 | MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
|
|---|
| 593 | MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
|
|---|
| 594 | fill1a.SetNameTab("PreCut");
|
|---|
| 595 | fill2a.SetNameTab("PostCut");
|
|---|
| 596 | fill3a.SetNameTab("VsSize");
|
|---|
| 597 | fill4a.SetNameTab("HilExt");
|
|---|
| 598 | fill5a.SetNameTab("HilSrc");
|
|---|
| 599 | fill6a.SetNameTab("ImgPar");
|
|---|
| 600 | fill7a.SetNameTab("NewPar");
|
|---|
| 601 |
|
|---|
| 602 | tlist2.AddToList(&read);
|
|---|
| 603 | tlist2.AddToList(&contsel);
|
|---|
| 604 | tlist2.AddToList(&calc);
|
|---|
| 605 | tlist2.AddToList(&hcalc1);
|
|---|
| 606 | tlist2.AddToList(&hcalc2);
|
|---|
| 607 | tlist2.AddToList(&fill1a);
|
|---|
| 608 | tlist2.AddToList(fCut0);
|
|---|
| 609 | tlist2.AddToList(fCut1);
|
|---|
| 610 | tlist2.AddToList(fCut2);
|
|---|
| 611 | tlist2.AddToList(fCut3);
|
|---|
| 612 | tlist2.AddToList(fEstimateEnergy);
|
|---|
| 613 | tlist2.AddToList(&fill3);
|
|---|
| 614 | tlist2.AddToList(&fill4);
|
|---|
| 615 | tlist2.AddToList(&fill2a);
|
|---|
| 616 | tlist2.AddToList(&fill3a);
|
|---|
| 617 | tlist2.AddToList(&fill4a);
|
|---|
| 618 | tlist2.AddToList(&fill5a);
|
|---|
| 619 | tlist2.AddToList(&fill6a);
|
|---|
| 620 | tlist2.AddToList(&fill7a);
|
|---|
| 621 |
|
|---|
| 622 | MEvtLoop loop2(fName);
|
|---|
| 623 | loop2.SetParList(&plist);
|
|---|
| 624 | loop2.SetDisplay(fDisplay);
|
|---|
| 625 | loop2.SetLogStream(fLog);
|
|---|
| 626 |
|
|---|
| 627 | if (!SetupEnv(loop2))
|
|---|
| 628 | return kFALSE;
|
|---|
| 629 |
|
|---|
| 630 | if (!loop2.Eventloop(fMaxEvents))
|
|---|
| 631 | {
|
|---|
| 632 | *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
|
|---|
| 633 | return kFALSE;
|
|---|
| 634 | }
|
|---|
| 635 |
|
|---|
| 636 | tlist2.PrintStatistics();
|
|---|
| 637 |
|
|---|
| 638 | if (!loop2.GetDisplay())
|
|---|
| 639 | {
|
|---|
| 640 | *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
|
|---|
| 641 | return kFALSE;
|
|---|
| 642 | }
|
|---|
| 643 |
|
|---|
| 644 | // -------------------------- Spectrum ----------------------------
|
|---|
| 645 |
|
|---|
| 646 | TH1D collarea(area.GetHEnergy());
|
|---|
| 647 | TH1D weights;
|
|---|
| 648 | hest.GetWeights(weights);
|
|---|
| 649 |
|
|---|
| 650 | cout << "Effective On time: " << ontime << "s" << endl;
|
|---|
| 651 |
|
|---|
| 652 | excess.SetDirectory(NULL);
|
|---|
| 653 | excess.SetBit(kCanDelete);
|
|---|
| 654 | excess.Scale(1./ontime);
|
|---|
| 655 | excess.Divide(&collarea);
|
|---|
| 656 | excess.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
|
|---|
| 657 | excess.SetYTitle("N/sm^{2}");
|
|---|
| 658 |
|
|---|
| 659 | TCanvas &c = fDisplay->AddTab("Spectrum");
|
|---|
| 660 | c.Divide(2,2);
|
|---|
| 661 | c.cd(1);
|
|---|
| 662 | gPad->SetBorderMode(0);
|
|---|
| 663 | gPad->SetLogx();
|
|---|
| 664 | gPad->SetLogy();
|
|---|
| 665 | gPad->SetGridx();
|
|---|
| 666 | gPad->SetGridy();
|
|---|
| 667 | collarea.DrawCopy();
|
|---|
| 668 |
|
|---|
| 669 | c.cd(2);
|
|---|
| 670 | gPad->SetBorderMode(0);
|
|---|
| 671 | gPad->SetLogx();
|
|---|
| 672 | gPad->SetLogy();
|
|---|
| 673 | gPad->SetGridx();
|
|---|
| 674 | gPad->SetGridy();
|
|---|
| 675 | excess.DrawCopy();
|
|---|
| 676 |
|
|---|
| 677 | c.cd(3);
|
|---|
| 678 | gPad->SetBorderMode(0);
|
|---|
| 679 | gPad->SetLogx();
|
|---|
| 680 | gPad->SetLogy();
|
|---|
| 681 | gPad->SetGridx();
|
|---|
| 682 | gPad->SetGridy();
|
|---|
| 683 | weights.DrawCopy();
|
|---|
| 684 |
|
|---|
| 685 | excess.Divide(&weights);
|
|---|
| 686 | //excess.Multiply(&weights);
|
|---|
| 687 | excess.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
|
|---|
| 688 |
|
|---|
| 689 | for (int i=0; i<excess.GetNbinsX(); i++)
|
|---|
| 690 | {
|
|---|
| 691 | excess.SetBinContent(i+1, excess.GetBinContent(i+1)/excess.GetBinWidth(i+1)*1000);
|
|---|
| 692 | excess.SetBinError(i+1, excess.GetBinError(i+1)/ excess.GetBinWidth(i+1)*1000);
|
|---|
| 693 | }
|
|---|
| 694 |
|
|---|
| 695 | excess.SetYTitle("N/TeVsm^{2}");
|
|---|
| 696 |
|
|---|
| 697 | c.cd(4);
|
|---|
| 698 | gPad->SetBorderMode(0);
|
|---|
| 699 | gPad->SetLogx();
|
|---|
| 700 | gPad->SetLogy();
|
|---|
| 701 | gPad->SetGridx();
|
|---|
| 702 | gPad->SetGridy();
|
|---|
| 703 | excess.DrawCopy();
|
|---|
| 704 |
|
|---|
| 705 | TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
|
|---|
| 706 | f.SetParameter(0, -2.87);
|
|---|
| 707 | f.SetParameter(1, 1.9e-6);
|
|---|
| 708 | f.SetLineColor(kGreen);
|
|---|
| 709 | excess.Fit(&f, "NI", "", 55, 2e4);
|
|---|
| 710 | f.DrawCopy("same");
|
|---|
| 711 |
|
|---|
| 712 | if (!fPathOut.IsNull())
|
|---|
| 713 | fDisplay->SaveAsRoot(fPathOut);
|
|---|
| 714 |
|
|---|
| 715 | /*
|
|---|
| 716 | TString str;
|
|---|
| 717 | str += "(1.68#pm0.15)10^{-7}";
|
|---|
| 718 | str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
|
|---|
| 719 | str += "\\frac{ph}{TeVm^{2}s}";
|
|---|
| 720 |
|
|---|
| 721 | TLatex tex;
|
|---|
| 722 | tex.DrawLatex(2e2, 7e-5, str);
|
|---|
| 723 | */
|
|---|
| 724 |
|
|---|
| 725 | return kTRUE;
|
|---|
| 726 | }
|
|---|