Changeset 7147


Ignore:
Timestamp:
06/13/05 09:41:57 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Cosy/tpoint/plot.C

    r7105 r7147  
    121121};
    122122
    123 const Int_t counts = 12-4;
     123const Int_t counts = 13-4;
    124124Description_t desc[counts] =
    125125{
    … …  
    139139    {"05051", "TPoints Residuals 5/2004-1" , "tpoint/tpoint0505-1.txt"},
    140140    // Mirror alignment has been fixed
    141     {"05052", "TPoints Residuals 5/2004-2" , "tpoint/tpoint0505-2.txt"}
     141    {"05052", "TPoints Residuals 5/2004-2" , "tpoint/tpoint0505-2.txt"},
     142    // Fixes to pointing model due to fixing a screw
     143    {"0506", "TPoints Residuals 6/2004" , "tpoint/tpoint0506.txt"}
    142144};
    143145

  • trunk/MagicSoft/Mars/Changelog

    r7146 r7147  
    2121
    2222                                                 -*-*- END OF LINE -*-*-
     23 2005/06/13 Thomas Bretz
     24
     25  * mcalib/MCalibrationHiLoCam.h:
     26     - added Print to //*MENU*
     27
     28   * mhflux/MHThetaSq.[h,cc]:
     29     - added resources for fNumBinsSignal and fNumBinsTotal
     30
     31   * mjobs/MJSpectrum.[h,cc]:
     32     - implemented weighting in theta, so we get better statistics
     33     - improved output
     34     - added plotting other spectras
     35
     36
     37
    2338 2005/06/10 Daniela Dorner
    2439

  • trunk/MagicSoft/Mars/mcalib/MCalibrationHiLoCam.h

    r5947 r7147  
    1919
    2020  // Prints
    21   void Print(Option_t *o="") const;
     21  void Print(Option_t *o="") const; //*MENU*
    2222
    2323  // Others

  • trunk/MagicSoft/Mars/mhflux/MHThetaSq.cc

    r7125 r7147  
    3434// For more detailes see MHAlpha.
    3535//
     36// Version 2:
     37// ---------
     38//  + UInt_t fNumBinsSignal;
     39//  + UInt_t fNumBinsTotal;
     40//
    3641//////////////////////////////////////////////////////////////////////////////
    3742#include "MHThetaSq.h"
    … …  
    5762//
    5863MHThetaSq::MHThetaSq(const char *name, const char *title)
    59     : MHAlpha(name, title), fThetaSq(0)
     64    : MHAlpha(name, title), fThetaSq(0), fNumBinsSignal(3), fNumBinsTotal(45)
    6065{
    6166    //
    … …  
    120125    // Calculate bining which fits alpha-cut
    121126    const Double_t intmax = fit->GetSignalIntegralMax();
    122     const UInt_t   nbins  = 75;
    123     const UInt_t   nsig   =  5;
     127    const UInt_t   nbins  = fNumBinsTotal;
     128    const UInt_t   nsig   = fNumBinsSignal;
    124129
    125130    MBinning binsa(nbins, 0, nbins*intmax/nsig);
    … …  
    201206   //     fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
    202207}
     208
     209Int_t MHThetaSq::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
     210{
     211    Int_t rc = MHAlpha::ReadEnv(env, prefix, print);
     212    if (rc==kERROR)
     213        return kERROR;
     214
     215    if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
     216    {
     217        SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
     218        rc = kTRUE;
     219    }
     220    if (IsEnvDefined(env, prefix, "NumBinsTotal", print))
     221    {
     222        SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal));
     223        rc = kTRUE;
     224    }
     225    return rc;
     226}

  • trunk/MagicSoft/Mars/mhflux/MHThetaSq.h

    r7125 r7147  
    1313    MParameterD  *fThetaSq; //!
    1414
     15    UInt_t fNumBinsSignal;
     16    UInt_t fNumBinsTotal;
     17
    1518    Bool_t      GetParameter(const MParList &pl);
    1619    Double_t    GetVal() const;
    … …  
    2326    void InitMapping(MHMatrix *mat, Int_t type=0);
    2427
     28    // MParContainer
     29    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
     30
    2531public:
    2632    MHThetaSq(const char *name=NULL, const char *title=NULL);
    2733
    28     ClassDef(MHThetaSq, 1) // Theta-Plot which is fitted online
     34    void SetNumBinsSignal(UInt_t n) { fNumBinsSignal=TMath::Max(n, 1U); }
     35    void SetNumBinsTotal(UInt_t n)  { fNumBinsTotal =TMath::Max(n, 1U); }
     36
     37    ClassDef(MHThetaSq, 2) // Theta-Plot which is fitted online
    2938};
    3039

  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

    r7099 r7147  
    3838#include <TFile.h>
    3939#include <TChain.h>
     40#include <TLatex.h>
    4041#include <TCanvas.h>
    4142#include <TObjArray.h>
    … …  
    8485MJSpectrum::MJSpectrum(const char *name, const char *title)
    8586    : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
    86     fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
     87    fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE),
     88    fNoThetaWeights(kFALSE)
    8789{
    8890    fName  = name  ? name  : "MJSpectrum";
    … …  
    336338    // Calculate the Probability
    337339    temp1.Divide(&temp2);
    338     temp1.Scale(1./temp1.GetMaximum());
     340    temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());
    339341
    340342    // Some cosmetics: Name, Axis, etc.
    … …  
    348350}
    349351
     352// --------------------------------------------------------------------------
     353//
     354// Display the final theta distribution.
     355//
    350356void MJSpectrum::DisplayResult(const TH2D &h2) const
    351357{
    … …  
    376382}
    377383
     384// --------------------------------------------------------------------------
     385//
     386// Fills the excess histogram (vs E-est) from the events stored in the
     387// ganymed result file and therefor estimates the energy.
     388//
     389// The resulting histogram excess-vs-energy ist copied into h2.
     390//
    378391Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
    379392{
    … …  
    531544}
    532545
     546// --------------------------------------------------------------------------
     547//
     548// Calculate the final spectrum from:
     549//  - collection area
     550//  - excess
     551//  - correction coefficients
     552//  - ontime
     553// and display it
     554//
    533555TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
    534556{
    … …  
    573595    weights.DrawCopy();
    574596
    575     //spectrum.Divide(&weights);
    576     spectrum.Multiply(&weights);
    577     spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
     597    //spectrum.Multiply(&weights);
     598    spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
     599    spectrum.SetBit(TH1::kNoStats);
    578600
    579601    for (int i=0; i<excess.GetNbinsX(); i++)
    580602    {
     603        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
     604        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
     605
    581606        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
    582607        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
    … …  
    589614    gPad->SetGridx();
    590615    gPad->SetGridy();
     616    spectrum.SetMinimum(1e-12);
    591617    spectrum.SetXTitle("E [GeV]");
    592618    spectrum.SetYTitle("N/TeVsm^{2}");
    593619    spectrum.DrawCopy();
    594620
    595     TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
     621    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
    596622    f.SetParameter(0, -2.87);
    597623    f.SetParameter(1, 1.9e-6);
    598624    f.SetLineColor(kGreen);
    599     spectrum.Fit(&f, "NIM", "", 55, 2e4);
     625    spectrum.Fit(&f, "NIM", "", 100, 4000);
    600626    f.DrawCopy("same");
    601627
    602     /*
    603      TString str;
    604      str += "(1.68#pm0.15)10^{-7}";
    605      str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
    606      str += "\\frac{ph}{TeVm^{2}s}";
    607 
    608      TLatex tex;
    609      tex.DrawLatex(2e2, 7e-5, str);
    610      */
     628    const Double_t p0 = f.GetParameter(0);
     629    const Double_t p1 = f.GetParameter(1);
     630
     631    const Double_t e0 = f.GetParError(0);
     632    const Double_t e1 = f.GetParError(1);
     633
     634    const Int_t    np  = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
     635    const Double_t exp = TMath::Power(10, np);
     636
     637    TString str;
     638    str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
     639    str += Form("(\\frac{E}{TeV})^{-%.2f#pm%.2f}", p0, e0);
     640    str += "\\frac{ph}{TeVm^{2}s}";
     641
     642    TLatex tex;
     643    tex.SetTextSize(0.045);
     644    tex.SetBit(TLatex::kTextNDC);
     645    tex.DrawLatex(0.45, 0.935, str);
     646
     647    str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
     648    tex.DrawLatex(0.70, 0.83, str);
    611649
    612650    TArrayD res(2);
    613651    res[0] = f.GetParameter(0);
    614652    res[1] = f.GetParameter(1);
     653
     654
     655    // Plot other spectra  from Whipple
     656    f.SetParameter(0, -2.45);
     657    f.SetParameter(1, 3.3e-7);
     658    f.SetRange(300, 8000);
     659    f.SetLineColor(kBlack);
     660    f.SetLineStyle(kDashed);
     661    f.DrawCopy("same");
     662
     663    // Plot other spectra  from Cangaroo
     664    f.SetParameter(0, -2.53);
     665    f.SetParameter(1, 2.0e-7);
     666    f.SetRange(7000, 50000);
     667    f.SetLineColor(kBlack);
     668    f.SetLineStyle(kDashed);
     669    f.DrawCopy("same");
     670
     671    // Plot other spectra  from Robert
     672    f.SetParameter(0, -2.59);
     673    f.SetParameter(1, 2.58e-7);
     674    f.SetRange(150, 1500);
     675    f.SetLineColor(kBlack);
     676    f.SetLineStyle(kDashed);
     677    f.DrawCopy("same");
     678
     679    // Plot other spectra  from HEGRA
     680    f.SetParameter(0, -2.61);
     681    f.SetParameter(1, 2.7e-7);
     682    f.SetRange(1000, 20000);
     683    f.SetLineColor(kBlack);
     684    f.SetLineStyle(kDashed);
     685    f.DrawCopy("same");
     686
    615687    return res;
    616688}
    617689
     690// --------------------------------------------------------------------------
     691//
     692// Scale some image parameter plots using the scale factor and plot them
     693// together with the corresponding MC histograms.
     694// Called from DisplaySize
     695//
    618696Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
    619697{
    … …  
    664742}
    665743
     744// --------------------------------------------------------------------------
     745//
     746// Take a lot of histograms and plot them together in one plot.
     747// Calls PlotSame
     748//
    666749Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
    667750{
    … …  
    698781    gPad->SetGridy();
    699782
    700     excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
     783    excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
    701784    excess = excess->DrawCopy("E2");
    702785    // Don't do this on the original object!
    … …  
    727810            histsel->KolmogorovTest(excess, "DX");
    728811            fLog->Separator("Chi^2 Test");
    729             histsel->Chi2Test(excess, "P");
     812            const Double_t p = histsel->Chi2Test(excess, "P");
     813
     814            TLatex tex;
     815            tex.SetBit(TLatex::kTextNDC);
     816            tex.DrawLatex(0.7, 0.93, Form("P(\\chi^{2})=%.0f", p*100));
    730817        }
    731818    }
    … …  
    814901    if (!GetThetaDistribution(temp1, temp2))
    815902        return kFALSE;
     903
     904    if (!fNoThetaWeights)
     905        weight.SetZdWeights(&temp1);
    816906
    817907    TH1D excess;
    … …  
    9301020
    9311021    tlist2.AddToList(&read);
    932     if (!fRawMc)
     1022    if (!fRawMc && fNoThetaWeights)
    9331023        tlist2.AddToList(&contsel);
    9341024    tlist2.AddToList(&calc);
Note: See TracChangeset for help on using the changeset viewer.