- Timestamp:
- 06/11/08 12:32:36 (16 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r8933 r8935 20 20 21 21 22 2008/06/11 Thomas Bretz 23 24 * mhflux/MHEnergyEst.[h,cc]: 25 - finally replaced fResolution by more correct histograms 26 - some code cleanup in projecting, profiling and drawing 27 - increased class version number by one 28 29 * mbase/MStatusDisplay.cc: 30 - remove the embedded canvas from the global list to prevent 31 global access to it 32 33 34 22 35 2008/06/10 Thomas Bretz 23 36 … … 41 54 - replaced some AddContainer by the new AddTree 42 55 - added Pyrometer information to in- and output, respectively 56 57 * datacenter/macros/fillstar.C: 58 - added new columns fAvgHumidity, fAvgCloudiness, fRmsCloudiness 59 and fAvgTempSky 60 61 * mhist/MHWeather.[h,cc]: 62 - removed the display of the solar radiation which was 63 never working 64 - added display of the pyrometer data to the display 65 - reorganized display 66 67 * mjobs/MJStar.cc: 68 - added filling of the weather data also from the pyrometer branch 43 69 44 70 -
trunk/MagicSoft/Mars/NEWS
r8933 r8935 61 61 * The code has been prepared for compilation with root 5.18/00d 62 62 63 * The MHEnergyEst histogram now shows the distribution of 64 (Eest-Emc)/Est and the distributions (Eest-Emc)/Eest vs. Eest 65 and (Eest-Emc)/Emc vs Emc. 66 63 67 ;merpp 64 68 … … 119 123 * The effective on-time calculation doesn't use events with only 120 124 sum-trigger anymore 125 126 * The data in the MHWeather tab has been reorganized. The never 127 working solar radiation has been removed and the data from 128 the pyrometer (cloudiness, air and sky temperature) is 129 displayed in addition. 121 130 122 131 ;ganymed/sponde -
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.