Changeset 8935 for trunk/MagicSoft/Mars/mhflux
- Timestamp:
- 06/11/08 12:32:36 (17 years ago)
- Location:
- trunk/MagicSoft/Mars/mhflux
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc
r8907 r8935 20 20 ! Author(s): Thomas Bretz 1/2005 <mailto:tbretz@astro.uni-wuerzburg.de> 21 21 ! 22 ! Copyright: MAGIC Software Development, 2000-200 522 ! Copyright: MAGIC Software Development, 2000-2008 23 23 ! 24 24 ! … … 31 31 // calculates the migration matrix E-est vs. E-true 32 32 // for different bins in Theta 33 // 34 // Class Version 2: 35 // - fHResolution 36 // + fHResolutionEst 37 // + fHResolutionMC 33 38 // 34 39 ////////////////////////////////////////////////////////////////////////////// … … 82 87 fHEnergy.SetZTitle("\\Theta [\\circ]"); 83 88 84 fHResolution.SetDirectory(NULL); 85 fHResolution.SetName("EnergyRes"); 86 fHResolution.SetTitle("Histogram in \\Delta E/E vs E_{est} and E_{mc}"); 87 fHResolution.SetXTitle("E_{est} [GeV]"); 88 fHResolution.SetYTitle("E_{mc} [GeV]"); 89 fHResolution.SetZTitle("E_{est}/E_{mc}-1"); 89 fHResolutionEst.SetDirectory(NULL); 90 fHResolutionEst.SetName("ResEnergyEst"); 91 fHResolutionEst.SetTitle("Histogram in \\Delta E/E_{est} vs E_{est}"); 92 fHResolutionEst.SetXTitle("E_{est} [GeV]"); 93 fHResolutionEst.SetYTitle("1-E_{mc}/E_{est}"); 94 95 fHResolutionMC.SetDirectory(NULL); 96 fHResolutionMC.SetName("ResEnergyMC"); 97 fHResolutionMC.SetTitle("Histogram in \\Delta E/E_{mc} vs E_{mc}"); 98 fHResolutionMC.SetXTitle("E_{mc} [GeV]"); 99 fHResolutionMC.SetYTitle("E_{est}/E_{mc}-1"); 90 100 91 101 fHImpact.SetDirectory(NULL); … … 105 115 SetBinning(&fHEnergy, &binse, &binse, &binst); 106 116 SetBinning(&fHImpact, &binsi, &binsr); 107 SetBinning(&fHResolution, &binse, &binse, &binsr); 117 118 SetBinning(&fHResolutionEst, &binse, &binsr); 119 SetBinning(&fHResolutionMC, &binse, &binsr); 108 120 109 121 // For some unknown reasons this must be called after 110 122 // the binning has been initialized at least once 111 123 fHEnergy.Sumw2(); 112 fHResolution.Sumw2();113 124 fHImpact.Sumw2(); 125 fHResolutionEst.Sumw2(); 126 fHResolutionMC.Sumw2(); 114 127 } 115 128 … … 154 167 SetBinning(&fHEnergy, &binse, &binse, &binst); 155 168 SetBinning(&fHImpact, &binsi, &binsr); 156 SetBinning(&fHResolution, &binse, &binse, &binsr); 169 170 SetBinning(&fHResolutionEst, &binse, &binsr); 171 SetBinning(&fHResolutionMC, &binse, &binsr); 157 172 158 173 fChisq = 0; … … 161 176 fHEnergy.Reset(); 162 177 fHImpact.Reset(); 163 fHResolution.Reset(); 178 179 fHResolutionEst.Reset(); 180 fHResolutionMC.Reset(); 164 181 165 182 return kTRUE; … … 176 193 const Double_t imp = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100; 177 194 const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg(); 178 const Double_t resE = (eest-etru)/etru; 195 196 const Double_t resEst = (eest-etru)/eest; 197 const Double_t resMC = (eest-etru)/etru; 179 198 180 199 fHEnergy.Fill(eest, etru, theta, w); 181 fHResolution.Fill(eest, etru, resE, w); 182 fHImpact.Fill(imp, resE, w); 200 fHImpact.Fill(imp, resEst, w); 201 202 fHResolutionEst.Fill(eest, resEst, w); 203 fHResolutionMC.Fill(etru, resMC, w); 183 204 184 205 // For the fit we use a different quantity … … 223 244 } 224 245 225 TH1D *h = (TH1D*)fHResolution .ProjectionZ("Resolution");246 TH1D *h = (TH1D*)fHResolutionEst.ProjectionY("Dummy", -1, -1, "s"); 226 247 h->Fit("gaus", "Q0", "", -1.0, 0.25); 227 248 … … 270 291 pad->cd(1); 271 292 272 TH1 D*hx=0;273 TH1 D*hy=0;293 TH1 *hx=0; 294 TH1 *hy=0; 274 295 275 296 if (pad->GetPad(1)) … … 277 298 pad->GetPad(1)->cd(1); 278 299 279 if ((hx=(TH1D*)gPad->FindObject("EnergyEst_ex"))) 280 { 281 TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ex"); 282 hx->Reset(); 283 hx->Add(h2); 284 delete h2; 285 } 286 287 if ((hy=(TH1D*)gPad->FindObject("EnergyEst_ey"))) 288 { 289 TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ey"); 290 hy->Reset(); 291 hy->Add(h2); 292 delete h2; 293 } 300 if (gPad->FindObject("EnergyEst_ex")) 301 hx = fHEnergy.Project3D("ex"); 302 303 if (gPad->FindObject("EnergyEst_ey")) 304 hy = fHEnergy.Project3D("ey"); 294 305 295 306 if (hx && hy) 296 307 { 308 hx->SetLineColor(kBlue); 309 hx->SetMarkerColor(kBlue); 297 310 hy->SetMaximum(); 298 311 hy->SetMaximum(TMath::Max(hx->GetMaximum(), hy->GetMaximum())*1.2); … … 305 318 pad->GetPad(1)->GetPad(2)->cd(1); 306 319 if ((hx=(TH1D*)gPad->FindObject("EnergyEst_ez"))) 307 { 308 TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ez"); 309 hx->Reset(); 310 hx->Add(h2); 311 delete h2; 312 } 320 fHEnergy.Project3D("ez"); 313 321 314 322 pad->GetPad(1)->GetPad(2)->cd(2); 315 hx = (TH1D*)fHResolution .ProjectionZ("Resolution");323 hx = (TH1D*)fHResolutionEst.ProjectionY("Resolution", -1, -1, "e"); 316 324 TPaveStats *stats = dynamic_cast<TPaveStats*>(hx->FindObject("stats")); 317 325 if (stats) … … 323 331 324 332 hx->Fit("gaus", "Q", "", -1.0, 0.25); 333 hx->GetFunction("gaus")->SetLineColor(kBlue); 334 hx->GetFunction("gaus")->SetLineWidth(2); 325 335 gPad=NULL; 326 336 gStyle->SetOptFit(101); … … 331 341 { 332 342 pad->GetPad(2)->cd(1); 333 UpdatePlot(fHEnergy, "yx", kTRUE); 343 if (gPad->FindObject("EnergyEst_yx")) 344 { 345 TH2D *hyx = static_cast<TH2D*>(fHEnergy.Project3D("yx")); 346 UpdateProf(*hyx, kTRUE); 347 } 334 348 335 349 TLine *l = (TLine*)gPad->FindObject("TLine"); … … 346 360 347 361 pad->GetPad(2)->cd(2); 348 UpdateP lot(fHResolution, "zy");362 UpdateProf(fHResolutionEst, kFALSE); 349 363 350 364 pad->GetPad(2)->cd(3); 351 UpdateP lot(fHResolution, "zx");365 UpdateProf(fHResolutionMC, kFALSE); 352 366 } 353 367 } 354 368 355 void MHEnergyEst::UpdatePlot(TH3 &h, const char *how, Bool_t logy) 356 { 357 TH2D *hyx=0; 358 if (!(hyx=(TH2D*)gPad->FindObject(MString::Format("%s_%s", h.GetName(), how)))) 369 void MHEnergyEst::UpdateProf(TH2 &h, Bool_t logy) 370 { 371 const TString pname = MString::Format("Prof%s", h.GetName()); 372 373 if (!gPad->FindObject(pname)) 359 374 return; 360 375 361 TH2D *h2 = (TH2D*)h.Project3D(MString::Format("dum_%s", how)); 362 hyx->Reset(); 363 hyx->Add(h2); 364 delete h2; 365 366 TH1D *hx = 0; 367 if ((hx=(TH1D*)gPad->FindObject(MString::Format("Prof%s", h.GetName())))) 368 { 369 hx = hyx->ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s"); 370 371 if (logy && hx->GetMaximum()>0) 372 gPad->SetLogy(); 373 } 374 } 375 376 TH1 *MHEnergyEst::MakePlot(TH3 &h, const char *how) 377 { 378 gPad->SetBorderMode(0); 376 TH1D *hx = h.ProfileX(pname, -1, -1, "s"); 377 hx->SetLineColor(kBlue); 378 hx->SetMarkerColor(kBlue); 379 380 if (logy && hx->GetMaximum()>0) 381 gPad->SetLogy(); 382 } 383 384 TH1 *MHEnergyEst::MakeProj(const char *how) 385 { 386 TH1 *p = fHEnergy.Project3D(how); 387 p->SetDirectory(NULL); 388 p->SetBit(kCanDelete); 389 p->SetBit(TH1::kNoStats); 390 p->SetMarkerStyle(kFullDotMedium); 391 p->SetLineColor(kBlue); 392 393 return p; 394 } 395 396 TH1 *MHEnergyEst::MakeProf(TH2 &h) 397 { 398 TH1 *p = h.ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s"); 399 p->SetDirectory(NULL); 400 p->SetBit(kCanDelete); 401 p->SetLineWidth(2); 402 p->SetLineColor(kBlue); 403 p->SetFillStyle(4000); 404 p->SetStats(kFALSE); 405 406 return p; 407 } 408 409 // -------------------------------------------------------------------------- 410 // 411 // Draw the histogram 412 // 413 void MHEnergyEst::Draw(Option_t *opt) 414 { 415 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); 416 417 // Do the projection before painting the histograms into 418 // the individual pads 419 AppendPad(""); 420 421 pad->SetBorderMode(0); 422 pad->Divide(2, 1, 1e-10, 1e-10); 423 424 TH1 *h; 425 426 // ---------------------------------------- 427 428 pad->cd(1); 429 gPad->SetBorderMode(0); 430 431 gPad->Divide(1, 2, 1e-10, 1e-10); 432 433 TVirtualPad *pad2 = gPad; 434 435 // ---------------------------------------- 436 437 pad2->cd(1); 438 gPad->SetBorderMode(0); 439 440 gPad->SetGridx(); 441 gPad->SetGridy(); 442 gPad->SetLogx(); 443 444 h = MakeProj("ey"); 445 h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)"); 446 h->SetXTitle("E [GeV]"); // E_mc 447 h->SetYTitle("Counts"); 448 h->Draw(); 449 450 h = MakeProj("ex"); 451 h->SetLineColor(kBlue); 452 h->SetMarkerColor(kBlue); 453 h->Draw("same"); 454 455 // ---------------------------------------- 456 457 pad2->cd(2); 458 gPad->SetBorderMode(0); 459 460 TVirtualPad *pad3 = gPad; 461 pad3->Divide(2, 1, 1e-10, 1e-10); 462 pad3->cd(1); 463 gPad->SetBorderMode(0); 464 gPad->SetGridx(); 465 gPad->SetGridy(); 466 467 h = MakeProj("ez"); 468 h->SetTitle("Zenith Angle Distribution"); 469 h->GetXaxis()->SetMoreLogLabels(); 470 h->GetXaxis()->SetNoExponent(); 471 h->Draw(); 472 473 // ---------------------------------------- 474 475 pad3->cd(2); 476 gPad->SetBorderMode(0); 477 gPad->SetGridx(); 478 gPad->SetGridy(); 479 480 h = fHResolutionEst.ProjectionY("_py"); 481 h->SetTitle("Distribution of \\Delta E/E_{est}"); 482 h->SetDirectory(NULL); 483 h->SetBit(kCanDelete); 484 h->GetXaxis()->SetRangeUser(-1.3, 1.3); 485 h->Draw(); 486 // ---------------------------------------- 487 488 pad->cd(2); 489 gPad->SetBorderMode(0); 490 491 gPad->Divide(1, 3, 1e-10, 1e-10); 492 pad2 = gPad; 493 494 // ---------------------------------------- 495 496 pad2->cd(1); 497 gPad->SetBorderMode(0); 498 gPad->SetLogy(); 379 499 gPad->SetLogx(); 380 500 gPad->SetGridx(); … … 384 504 //gROOT->GetListOfCleanups()->Add(gPad); // WHY? 385 505 386 TH2D *h2 = (TH2D*) h.Project3D(how);506 TH2D *h2 = (TH2D*)fHEnergy.Project3D("yx"); 387 507 h2->SetDirectory(NULL); 388 508 h2->SetBit(kCanDelete); 389 509 h2->SetFillColor(kBlue); 390 h2->SetLineColor(kRed); 391 392 TH1D *h1 = h2->ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s"); 393 h1->SetDirectory(NULL); 394 h1->SetBit(kCanDelete); 395 h1->SetLineWidth(2); 396 h1->SetLineColor(kBlue); 397 h1->SetFillStyle(4000); 398 h1->SetStats(kFALSE); 510 h2->SetTitle("Estimated energy E_{est} vs Monte Carlo energy E_{mc}"); 511 512 TH1 *h1 = MakeProf(*h2); 399 513 400 514 h2->Draw(""); 401 h1->Draw("E0 hist C same"); 402 // h1->Draw("Chistsame"); 403 404 return h2; 405 } 406 407 408 // -------------------------------------------------------------------------- 409 // 410 // Draw the histogram 411 // 412 void MHEnergyEst::Draw(Option_t *opt) 413 { 414 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); 415 416 // Do the projection before painting the histograms into 417 // the individual pads 418 AppendPad(""); 419 420 pad->SetBorderMode(0); 421 pad->Divide(2, 1, 1e-10, 1e-10); 422 423 TH1 *h; 424 425 pad->cd(1); 426 gPad->SetBorderMode(0); 427 428 gPad->Divide(1, 2, 1e-10, 1e-10); 429 430 TVirtualPad *pad2 = gPad; 431 432 pad2->cd(1); 433 gPad->SetBorderMode(0); 434 435 gPad->SetGridx(); 436 gPad->SetGridy(); 437 gPad->SetLogx(); 438 h = (TH1D*)fHEnergy.Project3D("ey"); 439 h->SetBit(TH1::kNoStats); 440 h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)"); 441 h->SetXTitle("E [GeV]"); // E_mc 442 h->SetYTitle("Counts"); 443 h->SetBit(kCanDelete); 444 h->SetDirectory(NULL); 445 h->SetMarkerStyle(kFullDotMedium); 446 h->Draw(); 447 448 h = (TH1D*)fHEnergy.Project3D("ex"); 449 h->SetBit(TH1::kNoTitle|TH1::kNoStats); 450 h->SetXTitle("E [GeV]"); // E_est 451 h->SetYTitle("Counts"); 452 h->SetBit(kCanDelete); 453 h->SetDirectory(NULL); 454 h->SetMarkerStyle(kFullDotMedium); 455 h->SetLineColor(kBlue); 456 h->SetMarkerColor(kBlue); 457 h->Draw("same"); 458 459 // FIXME: LEGEND 460 461 pad2->cd(2); 462 gPad->SetBorderMode(0); 463 464 TVirtualPad *pad3 = gPad; 465 pad3->Divide(2, 1, 1e-10, 1e-10); 466 pad3->cd(1); 467 gPad->SetBorderMode(0); 468 gPad->SetGridx(); 469 gPad->SetGridy(); 470 h = fHEnergy.Project3D("ez"); 471 h->SetTitle("Zenith Angle Distribution"); 472 h->SetBit(TH1::kNoStats); 473 h->SetDirectory(NULL); 474 h->SetXTitle("\\Theta [\\circ]"); 475 h->SetBit(kCanDelete); 476 h->Draw(); 477 478 pad3->cd(2); 479 gPad->SetBorderMode(0); 480 gPad->SetGridx(); 481 gPad->SetGridy(); 482 /* 483 h = fHImpact.ProjectionX("Impact", -1, -1, "e"); 484 h->SetBit(TH1::kNoStats); 485 h->SetTitle("Distribution of Impact"); 486 h->SetDirectory(NULL); 487 h->SetXTitle("Impact [m]"); 488 h->SetBit(kCanDelete); 489 h->Draw();*/ 490 // ---------------------------------------- 491 h = fHResolution.ProjectionZ("Resolution"); 492 //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); 493 h->SetYTitle("Counts"); 494 h->SetTitle("Distribution of \\Delta E/E"); 495 h->SetDirectory(NULL); 496 h->SetBit(kCanDelete); 497 h->GetXaxis()->SetRangeUser(-1.3, 1.3); 498 h->Draw(""); 499 //h->Fit("gaus"); 500 // ---------------------------------------- 501 502 pad->cd(2); 503 gPad->SetBorderMode(0); 504 505 gPad->Divide(1, 3, 1e-10, 1e-10); 506 pad2 = gPad; 507 508 pad2->cd(1); 509 gPad->SetLogy(); 510 h = MakePlot(fHEnergy, "xy"); 511 h->SetXTitle("E_{mc} [GeV]"); 512 h->SetYTitle("E_{est} [GeV]"); 513 h->SetTitle("Estimated energy E_{est} vs Monte Carlo energy E_{mc}"); 515 h1->Draw("E0 same"); 514 516 515 517 TLine line; … … 520 522 line.SetLineStyle(kDashed); 521 523 524 // ---------------------------------------- 525 522 526 pad2->cd(2); 523 h = MakePlot(fHResolution, "zy"); 524 h->SetXTitle("E_{mc} [GeV]"); 525 h->SetYTitle("(E_{est}/E_{mc}-1"); 526 h->SetTitle("Energy resolution \\Delta E/E vs Monte Carlo Energy E_{mc}"); 527 h->SetMinimum(-1.3); 528 h->SetMaximum(1.3); 529 530 line.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0); 527 gPad->SetBorderMode(0); 528 gPad->SetLogx(); 529 gPad->SetGridx(); 530 gPad->SetGridy(); 531 fHResolutionEst.Draw(); 532 MakeProf(fHResolutionEst)->Draw("E0 same"); 533 534 fHResolutionEst.GetXaxis()->SetMoreLogLabels(); 535 fHResolutionEst.GetXaxis()->SetNoExponent(); 536 537 line.DrawLine(fHResolutionEst.GetXaxis()->GetXmin(), 0, fHResolutionEst.GetXaxis()->GetXmax(), 0); 538 539 // ---------------------------------------- 531 540 532 541 pad2->cd(3); 533 h = MakePlot(fHResolution, "zx"); 534 h->SetXTitle("E_{est} [GeV]"); 535 h->SetYTitle("(E_{est}/E_{mc}-1"); 536 h->SetTitle("Energy resolution \\Delta E/E vs estimated Energy E_{est}"); 537 h->SetMinimum(-1.3); 538 h->SetMaximum(1.3); 539 540 line.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0); 542 gPad->SetBorderMode(0); 543 gPad->SetLogx(); 544 gPad->SetGridx(); 545 gPad->SetGridy(); 546 fHResolutionMC.Draw(); 547 MakeProf(fHResolutionMC)->Draw("E0 same"); 548 549 fHResolutionMC.GetXaxis()->SetMoreLogLabels(); 550 fHResolutionMC.GetXaxis()->SetNoExponent(); 551 552 line.DrawLine(fHResolutionMC.GetXaxis()->GetXmin(), 0, fHResolutionMC.GetXaxis()->GetXmax(), 0); 541 553 } 542 554 -
trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h
r7692 r8935 33 33 34 34 TH3D fHEnergy; 35 TH3D fHResolution; 35 TH2D fHResolutionEst; 36 TH2D fHResolutionMC; 36 37 TH2D fHImpact; 37 38 … … 39 40 Double_t fBias; 40 41 41 TH1 *MakePlot(TH3 &h, const char *how); 42 void UpdatePlot(TH3 &h, const char *how, Bool_t logy=kFALSE); 42 TH1 *MakeProj(const char *how); 43 TH1 *MakeProf(TH2 &h); 44 void UpdateProf(TH2 &h, Bool_t logy); 43 45 44 46 Double_t GetVal(Int_t i) const; … … 62 64 void Print(Option_t *o="") const; 63 65 64 //ClassDef(MHEnergyEst, 2) // 65 ClassDef(MHEnergyEst, 1) // Histogram for the result of the energy reconstruction 66 ClassDef(MHEnergyEst, 2) // Histogram for the result of the energy reconstruction 66 67 }; 67 68
Note:
See TracChangeset
for help on using the changeset viewer.