Commit 696c32c22a98e6dabcf1947d96b2b8345f7fc7b2
1 parent
6b6b5700
Clean up code
Showing
1 changed file
with
0 additions
and
119 deletions
scripts/cta_kb_root2caldb.py
| ... | ... | @@ -397,9 +397,6 @@ class caldb(): |
| 397 | 397 | # Write header |
| 398 | 398 | for file in files: |
| 399 | 399 | file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") |
| 400 | - file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") | |
| 401 | - file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") | |
| 402 | - file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") | |
| 403 | 400 | |
| 404 | 401 | # Get histograms. Konrad's root files do not have an EffectiveArea80 |
| 405 | 402 | # histogram, but his effective area is for 80% containment radius |
| ... | ... | @@ -502,122 +499,6 @@ class caldb(): |
| 502 | 499 | return |
| 503 | 500 | |
| 504 | 501 | |
| 505 | -# ============== # | |
| 506 | -# Generate files # | |
| 507 | -# ============== # | |
| 508 | -def make_files(filename, prfname, bkgname): | |
| 509 | - """ | |
| 510 | - Generate files. | |
| 511 | - | |
| 512 | - Parameters: | |
| 513 | - filename - Performance file name | |
| 514 | - prfname - Performance file name | |
| 515 | - bkgname - File function filename | |
| 516 | - """ | |
| 517 | - # Open performance file and file function | |
| 518 | - file = TFile(filename) | |
| 519 | - fprf = open(prfname, "w") | |
| 520 | - fbkg = open(bkgname, "w") | |
| 521 | - | |
| 522 | - # Write header | |
| 523 | - fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") | |
| 524 | - #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") | |
| 525 | - | |
| 526 | - # Get histograms. Konrad's root files do not have an EffectiveArea80 | |
| 527 | - # histogram, but his effective area is for 80% containment radius | |
| 528 | - # anyways, so we simply read this histogram into aeff80. | |
| 529 | - sens = TH1F() | |
| 530 | - bgrate = TH1F() | |
| 531 | - bgratesqdeg = TH1F() | |
| 532 | - aeff = TH1F() | |
| 533 | - aeff80 = TH1F() | |
| 534 | - angres = TH1F() | |
| 535 | - angres80 = TH1F() | |
| 536 | - eres = TH1F() | |
| 537 | - file.GetObject("DiffSens", sens) | |
| 538 | - file.GetObject("BGRate", bgrate) | |
| 539 | - file.GetObject("BGRatePerSqDeg", bgratesqdeg) | |
| 540 | - file.GetObject("EffectiveArea", aeff) | |
| 541 | - try: | |
| 542 | - file.GetObject("EffectiveArea80", aeff80) | |
| 543 | - except: | |
| 544 | - file.GetObject("EffectiveArea", aeff80) | |
| 545 | - file.GetObject("AngRes", angres) | |
| 546 | - file.GetObject("AngRes80", angres80) | |
| 547 | - file.GetObject("ERes", eres) | |
| 548 | - | |
| 549 | - # Setup vectors | |
| 550 | - energies = bgrate.GetXaxis() | |
| 551 | - nbins_eng = energies.GetNbins() | |
| 552 | - for i in range(nbins_eng): | |
| 553 | - | |
| 554 | - # Get logE and background rate per squared degree | |
| 555 | - logE = energies.GetBinCenter(i+1) | |
| 556 | - bgd = bgrate.GetBinContent(i+1,1) | |
| 557 | - bgdsq = bgratesqdeg.GetBinContent(i+1,1) | |
| 558 | - area = aeff.GetBinContent(i+1,1) | |
| 559 | - area80 = aeff80.GetBinContent(i+1,1) | |
| 560 | - r68 = angres.GetBinContent(i+1,1) | |
| 561 | - r80 = angres80.GetBinContent(i+1,1) | |
| 562 | - dE = eres.GetBinContent(i+1,1) | |
| 563 | - diffs = sens.GetBinContent(i+1,1) | |
| 564 | - | |
| 565 | - # Skip any NaNs | |
| 566 | - if math.isnan(area80): | |
| 567 | - continue | |
| 568 | - | |
| 569 | - # Compute energy (in MeV) | |
| 570 | - energy = math.pow(10.0, logE)*1.0e6 | |
| 571 | - emin = math.pow(10.0, logE-0.1)*1.0e6 | |
| 572 | - emax = math.pow(10.0, logE+0.1)*1.0e6 | |
| 573 | - ewidth = emax-emin | |
| 574 | - | |
| 575 | - # Convert into background rate per steradian and MeV | |
| 576 | - bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth | |
| 577 | - | |
| 578 | - # Compute control background rate (this only works for KB files | |
| 579 | - # as they are for 80% containment radius). But we don't use this | |
| 580 | - # anyways, it's just for control on display ... | |
| 581 | - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) | |
| 582 | - if omega > 0.0: | |
| 583 | - bkg_rate2 = bgd / omega / ewidth | |
| 584 | - else: | |
| 585 | - bkg_rate2 = 0.0 | |
| 586 | - | |
| 587 | - # Compute full effective area by dividing area80/0.80 | |
| 588 | - area = area80/0.80 | |
| 589 | - | |
| 590 | - # Write results in file | |
| 591 | - line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ | |
| 592 | - (logE, area, r68, r80, dE, bgd, diffs) | |
| 593 | - fprf.write(line+"\n") | |
| 594 | - fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") | |
| 595 | - | |
| 596 | - # Show performance file | |
| 597 | - #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \ | |
| 598 | - # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2) | |
| 599 | - | |
| 600 | - # Write trailer | |
| 601 | - fprf.write("---------------------------------------------\n") | |
| 602 | - fprf.write("Notes\n") | |
| 603 | - fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n") | |
| 604 | - fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") | |
| 605 | - fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") | |
| 606 | - fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") | |
| 607 | - fprf.write(" 5) Fractional Energy Resolution (rms)\n") | |
| 608 | - fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") | |
| 609 | - fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") | |
| 610 | - fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") | |
| 611 | - fprf.write(" systematics and statistics and at least 10 photons.\n") | |
| 612 | - | |
| 613 | - # Close files | |
| 614 | - fprf.close() | |
| 615 | - fbkg.close() | |
| 616 | - | |
| 617 | - # Return | |
| 618 | - return | |
| 619 | - | |
| 620 | - | |
| 621 | 502 | #==========================# |
| 622 | 503 | # Main routine entry point # |
| 623 | 504 | #==========================# | ... | ... |