| 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): Wolfgang Wittek 4/2002 <mailto:wittek@mppmu.mpg.de>
|
|---|
| 19 | ! Abelardo Moralejo 5/2003 <mailto:moralejo@pd.infn.it>
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2002
|
|---|
| 22 | !
|
|---|
| 23 | !
|
|---|
| 24 | \* ======================================================================== */
|
|---|
| 25 |
|
|---|
| 26 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 27 | // //
|
|---|
| 28 | // MHMcEnergyMigration //
|
|---|
| 29 | // //
|
|---|
| 30 | // calculates the migration matrix E-est vs. E-true //
|
|---|
| 31 | // for different bins in Theta //
|
|---|
| 32 | // //
|
|---|
| 33 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 34 |
|
|---|
| 35 | #include "MHMcEnergyMigration.h"
|
|---|
| 36 |
|
|---|
| 37 | #include <TCanvas.h>
|
|---|
| 38 |
|
|---|
| 39 | #include "MMcEvt.hxx"
|
|---|
| 40 |
|
|---|
| 41 | #include "MEnergyEst.h"
|
|---|
| 42 | #include "MBinning.h"
|
|---|
| 43 | #include "MHillasSrc.h"
|
|---|
| 44 | #include "MParList.h"
|
|---|
| 45 |
|
|---|
| 46 | #include "MLog.h"
|
|---|
| 47 | #include "MLogManip.h"
|
|---|
| 48 | #include <TProfile.h>
|
|---|
| 49 |
|
|---|
| 50 | ClassImp(MHMcEnergyMigration);
|
|---|
| 51 |
|
|---|
| 52 |
|
|---|
| 53 | // --------------------------------------------------------------------------
|
|---|
| 54 | //
|
|---|
| 55 | // Default Constructor. It sets name and title of the histogram.
|
|---|
| 56 | //
|
|---|
| 57 | MHMcEnergyMigration::MHMcEnergyMigration(const char *name, const char *title)
|
|---|
| 58 | : fHist(), fHist2(), fHistImp()
|
|---|
| 59 | {
|
|---|
| 60 | //
|
|---|
| 61 | // set the name and title of this object
|
|---|
| 62 | //
|
|---|
| 63 | fName = name ? name : "MHMcEnergyMigration";
|
|---|
| 64 | fTitle = title ? title : "3-D histogram E-true E-est Theta";
|
|---|
| 65 |
|
|---|
| 66 | fHist.SetDirectory(NULL);
|
|---|
| 67 | fHist.SetTitle("3D-plot E_{TRUE} E_{EST} \\Theta");
|
|---|
| 68 | fHist.SetXTitle("E_{TRUE} [GeV]");
|
|---|
| 69 | fHist.SetYTitle("E_{EST} [GeV]");
|
|---|
| 70 | fHist.SetZTitle("\\Theta [\\circ]");
|
|---|
| 71 |
|
|---|
| 72 | fHist2.SetDirectory(NULL);
|
|---|
| 73 | fHist2.SetTitle("3D-plot \\Delta E / E vs E_{TRUE} E_{EST}");
|
|---|
| 74 | fHist2.SetXTitle("E_{TRUE} [GeV]");
|
|---|
| 75 | fHist2.SetYTitle("E_{EST} [GeV]");
|
|---|
| 76 | fHist2.SetZTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
|
|---|
| 77 |
|
|---|
| 78 | fHistImp.SetDirectory(NULL);
|
|---|
| 79 | fHistImp.SetTitle("\\Delta E / E vs Impact parameter");
|
|---|
| 80 | fHistImp.SetXTitle("Impact parameter (m)");
|
|---|
| 81 | fHistImp.SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
|
|---|
| 82 |
|
|---|
| 83 | }
|
|---|
| 84 |
|
|---|
| 85 | // --------------------------------------------------------------------------
|
|---|
| 86 | //
|
|---|
| 87 | // Set the binnings and prepare the filling of the histograms
|
|---|
| 88 | //
|
|---|
| 89 | Bool_t MHMcEnergyMigration::SetupFill(const MParList *plist)
|
|---|
| 90 | {
|
|---|
| 91 | fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst");
|
|---|
| 92 | if (!fEnergy)
|
|---|
| 93 | {
|
|---|
| 94 | *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
|
|---|
| 95 | return kFALSE;
|
|---|
| 96 | }
|
|---|
| 97 |
|
|---|
| 98 | fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
|
|---|
| 99 | if (!fMcEvt)
|
|---|
| 100 | {
|
|---|
| 101 | *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
|
|---|
| 102 | return kFALSE;
|
|---|
| 103 | }
|
|---|
| 104 |
|
|---|
| 105 | const MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
|
|---|
| 106 | const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
|
|---|
| 107 | if (!binsenergy || !binstheta)
|
|---|
| 108 | {
|
|---|
| 109 | *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
|
|---|
| 110 | return kFALSE;
|
|---|
| 111 | }
|
|---|
| 112 |
|
|---|
| 113 | SetBinning(&fHist, binsenergy, binsenergy, binstheta);
|
|---|
| 114 | fHist.Sumw2();
|
|---|
| 115 |
|
|---|
| 116 | const MBinning* binsde = (MBinning*)plist->FindObject("BinningDE");
|
|---|
| 117 | const MBinning* binsimpact = (MBinning*)plist->FindObject("BinningImpact");
|
|---|
| 118 | SetBinning(&fHistImp, binsimpact, binsde);
|
|---|
| 119 |
|
|---|
| 120 | fHistImp.Sumw2();
|
|---|
| 121 |
|
|---|
| 122 | SetBinning(&fHist2, binsenergy, binsenergy, binsde);
|
|---|
| 123 | fHist2.Sumw2();
|
|---|
| 124 |
|
|---|
| 125 | return kTRUE;
|
|---|
| 126 | }
|
|---|
| 127 |
|
|---|
| 128 | // --------------------------------------------------------------------------
|
|---|
| 129 | //
|
|---|
| 130 | // Draw the histogram
|
|---|
| 131 | //
|
|---|
| 132 | void MHMcEnergyMigration::Draw(Option_t *opt)
|
|---|
| 133 | {
|
|---|
| 134 | TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
|
|---|
| 135 | pad->SetBorderMode(0);
|
|---|
| 136 |
|
|---|
| 137 | AppendPad("");
|
|---|
| 138 |
|
|---|
| 139 | pad->Divide(2,2);
|
|---|
| 140 |
|
|---|
| 141 | TH1 *h;
|
|---|
| 142 |
|
|---|
| 143 | pad->cd(1);
|
|---|
| 144 | gPad->SetBorderMode(0);
|
|---|
| 145 | gPad->SetLogx();
|
|---|
| 146 | h = fHist.Project3D("ex_pro");
|
|---|
| 147 | h->SetTitle("Distribution of E_{TRUE}");
|
|---|
| 148 | h->SetXTitle("E_{TRUE} [GeV]");
|
|---|
| 149 | h->SetBit(kCanDelete);
|
|---|
| 150 | h->Draw(opt);
|
|---|
| 151 | h->SetDirectory(NULL);
|
|---|
| 152 |
|
|---|
| 153 | pad->cd(2);
|
|---|
| 154 | gPad->SetBorderMode(0);
|
|---|
| 155 | gPad->SetLogx();
|
|---|
| 156 | h = fHist.Project3D("ey_pro");
|
|---|
| 157 | h->SetTitle("Distribution of E_{EST}");
|
|---|
| 158 | h->SetXTitle("E_{est} [GeV]");
|
|---|
| 159 | h->SetBit(kCanDelete);
|
|---|
| 160 | h->Draw(opt);
|
|---|
| 161 | Double_t minEest = h->GetXaxis()->GetXmin();
|
|---|
| 162 | h->SetDirectory(NULL);
|
|---|
| 163 |
|
|---|
| 164 | pad->cd(3);
|
|---|
| 165 | gPad->SetBorderMode(0);
|
|---|
| 166 | h = fHist.Project3D("z_pro");
|
|---|
| 167 | h->SetTitle("Distribution of \\Theta");
|
|---|
| 168 | h->SetXTitle("\\Theta [\\circ]");
|
|---|
| 169 | h->SetBit(kCanDelete);
|
|---|
| 170 | h->SetLineWidth(2);
|
|---|
| 171 | h->Draw(opt);
|
|---|
| 172 | h->SetDirectory(NULL);
|
|---|
| 173 |
|
|---|
| 174 | pad->cd(4);
|
|---|
| 175 | gPad->SetBorderMode(0);
|
|---|
| 176 | TH1D* hpx;
|
|---|
| 177 | hpx = fHistImp.ProjectionX("_px", 1, fHistImp.GetNbinsY(),"e");
|
|---|
| 178 | hpx->SetTitle("Distribution of Impact parameter");
|
|---|
| 179 | hpx->SetXTitle("Impact parameter (m)");
|
|---|
| 180 | hpx->SetBit(kCanDelete);
|
|---|
| 181 | hpx->Draw(opt);
|
|---|
| 182 | hpx->SetDirectory(NULL);
|
|---|
| 183 | fHistImp.SetDirectory(NULL);
|
|---|
| 184 |
|
|---|
| 185 | pad->Modified();
|
|---|
| 186 | pad->Update();
|
|---|
| 187 |
|
|---|
| 188 | TVirtualPad *pad2 = MakeDefCanvas((TString)GetName()+"2");
|
|---|
| 189 | pad2->SetBorderMode(0);
|
|---|
| 190 |
|
|---|
| 191 | AppendPad("");
|
|---|
| 192 |
|
|---|
| 193 | pad2->Divide(2,2);
|
|---|
| 194 |
|
|---|
| 195 | TH2D *h2;
|
|---|
| 196 |
|
|---|
| 197 | pad2->cd(1);
|
|---|
| 198 | gPad->SetBorderMode(0);
|
|---|
| 199 | gPad->SetLogx();
|
|---|
| 200 | gPad->SetLogy();
|
|---|
| 201 | h2 = (TH2D*) fHist.Project3D("xy");
|
|---|
| 202 |
|
|---|
| 203 | TProfile* h2pfx;
|
|---|
| 204 | h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S");
|
|---|
| 205 | h2pfx->SetXTitle("E_{TRUE} (GeV)");
|
|---|
| 206 | h2pfx->SetYTitle("E_{EST} (GeV)");
|
|---|
| 207 | h2pfx->SetTitle("E_{EST} vs E_{TRUE}");
|
|---|
| 208 | h2pfx->SetBit(kCanDelete);
|
|---|
| 209 | h2pfx->SetFillColor(41);
|
|---|
| 210 | h2pfx->SetMinimum(minEest);
|
|---|
| 211 | h2pfx->SetStats(kFALSE);
|
|---|
| 212 | h2pfx->Draw("E3");
|
|---|
| 213 |
|
|---|
| 214 | h2->SetBit(kCanDelete);
|
|---|
| 215 | h2->SetFillColor(1);
|
|---|
| 216 | h2->Draw("box,same");
|
|---|
| 217 | h2->SetDirectory(NULL);
|
|---|
| 218 |
|
|---|
| 219 | h2pfx->SetLineColor(2);
|
|---|
| 220 | h2pfx->SetLineWidth(2);
|
|---|
| 221 | h2pfx->Draw("L,histo,same");
|
|---|
| 222 | h2pfx->SetDirectory(NULL);
|
|---|
| 223 |
|
|---|
| 224 | pad2->cd(2);
|
|---|
| 225 | gPad->SetBorderMode(0);
|
|---|
| 226 | gPad->SetLogx();
|
|---|
| 227 | gPad->SetLeftMargin(0.15);
|
|---|
| 228 | h2 = (TH2D*) fHist2.Project3D("zx");
|
|---|
| 229 | h2->SetBit(kCanDelete);
|
|---|
| 230 | h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S");
|
|---|
| 231 | h2pfx->SetTitle("\\Delta E / E vs E_{TRUE}");
|
|---|
| 232 | h2pfx->SetXTitle("E_{TRUE} (GeV)");
|
|---|
| 233 | h2pfx->SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
|
|---|
| 234 | h2pfx->SetBit(kCanDelete);
|
|---|
| 235 | h2pfx->SetFillColor(41);
|
|---|
| 236 | h2pfx->GetYaxis()->SetTitleOffset(1.4);
|
|---|
| 237 | h2pfx->SetStats(kFALSE);
|
|---|
| 238 | h2pfx->Draw("E3");
|
|---|
| 239 | h2->SetFillColor(1);
|
|---|
| 240 | h2->Draw("same,box");
|
|---|
| 241 | h2->SetDirectory(NULL);
|
|---|
| 242 | h2pfx->SetLineColor(2);
|
|---|
| 243 | h2pfx->SetLineWidth(2);
|
|---|
| 244 | h2pfx->Draw("L,histo,same");
|
|---|
| 245 | h2pfx->SetDirectory(NULL);
|
|---|
| 246 |
|
|---|
| 247 | pad2->cd(3);
|
|---|
| 248 | gPad->SetBorderMode(0);
|
|---|
| 249 | gPad->SetLogx();
|
|---|
| 250 | gPad->SetLeftMargin(0.15);
|
|---|
| 251 | h2 = (TH2D*) fHist2.Project3D("zy");
|
|---|
| 252 | h2->SetBit(kCanDelete);
|
|---|
| 253 | h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S");
|
|---|
| 254 | h2pfx->SetTitle("\\Delta E / E vs E_{EST}");
|
|---|
| 255 | h2pfx->SetXTitle("E_{EST} (GeV)");
|
|---|
| 256 | h2pfx->SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
|
|---|
| 257 | h2pfx->SetBit(kCanDelete);
|
|---|
| 258 | h2pfx->SetFillColor(41);
|
|---|
| 259 | h2pfx->GetYaxis()->SetTitleOffset(1.4);
|
|---|
| 260 | h2pfx->SetStats(kFALSE);
|
|---|
| 261 | h2pfx->SetMinimum(-1.);
|
|---|
| 262 | h2pfx->SetMaximum(1.);
|
|---|
| 263 | h2pfx->Draw("E3");
|
|---|
| 264 |
|
|---|
| 265 | h2->SetFillColor(1);
|
|---|
| 266 | h2->Draw("same,box");
|
|---|
| 267 | h2->SetDirectory(NULL);
|
|---|
| 268 | h2pfx->SetLineColor(2);
|
|---|
| 269 | h2pfx->SetLineWidth(2);
|
|---|
| 270 | h2pfx->Draw("L,histo,same");
|
|---|
| 271 | h2pfx->SetDirectory(NULL);
|
|---|
| 272 |
|
|---|
| 273 | pad2->cd(4);
|
|---|
| 274 | gPad->SetBorderMode(0);
|
|---|
| 275 | h = fHist2.ProjectionZ("_pz",1,fHist2.GetNbinsX(),1,fHist2.GetNbinsY(),"e");
|
|---|
| 276 | h->SetBit(kCanDelete);
|
|---|
| 277 | h->Draw();
|
|---|
| 278 | h->SetDirectory(NULL);
|
|---|
| 279 |
|
|---|
| 280 | pad2->Modified();
|
|---|
| 281 | pad2->Update();
|
|---|
| 282 |
|
|---|
| 283 | fHist.SetDirectory(NULL);
|
|---|
| 284 | fHist2.SetDirectory(NULL);
|
|---|
| 285 |
|
|---|
| 286 | }
|
|---|
| 287 |
|
|---|
| 288 | // --------------------------------------------------------------------------
|
|---|
| 289 | //
|
|---|
| 290 | // Fill the histogram
|
|---|
| 291 | //
|
|---|
| 292 | Bool_t MHMcEnergyMigration::Fill(const MParContainer *par, const Stat_t w)
|
|---|
| 293 | {
|
|---|
| 294 | // get E-true from fMcEvt and E-est from fEnergy
|
|---|
| 295 |
|
|---|
| 296 | fHist.Fill(fMcEvt->GetEnergy(), fEnergy->GetEnergy(), fMcEvt->GetTelescopeTheta()*kRad2Deg, w);
|
|---|
| 297 |
|
|---|
| 298 | fHist2.Fill(fMcEvt->GetEnergy(), fEnergy->GetEnergy(), (fEnergy->GetEnergy()-fMcEvt->GetEnergy())/fMcEvt->GetEnergy(), w);
|
|---|
| 299 |
|
|---|
| 300 | fHistImp.Fill(fMcEvt->GetImpact()/100., (fEnergy->GetEnergy()-fMcEvt->GetEnergy())/fMcEvt->GetEnergy(), w);
|
|---|
| 301 |
|
|---|
| 302 | return kTRUE;
|
|---|
| 303 | }
|
|---|
