CTAClassicalAnalysis.py

Martin Pierrick, 06/23/2013 04:11 PM

Download (3.53 KB)

 
1 
def ClassicalAnalysis():
2 
    """
3 
    Typical script for classical analysis of TeV data
4 
    (virtual at this stage, only to help organizing the implementation)
5
6 
    Parameters
7 
    - ....
8 
    - ....
9
10 
    Notes
11 
    - ...
12 
    - ...
13 
    """
14
15
16 
    # Load/simulate data
17 
    ####################
18
19 
    # Load data if they already exist somewhere...
20 
    if (os.path.exists('obs.xml')):
21 
            obs = GObservations('obs.xml')
22 
    #... else simulate them
23 
    else:
24 
            obs = ...
25
26 
    # We can inspect data
27 
    for run in obs:
28 
            run.print()
29
30 
    # Set analysis parameters
31 
    #########################
32
33 
    # Define energy binning
34 
    ebds = GEbounds(10, GEnergy(src_emin, "TeV"), GEnergy(src_emax, "TeV"))
35
36 
    # Define ON region
37 
    on_centre=GSkyDir(...)
38 
    on_radius=0.1
39 
    on_reg=GSkyRegionCircle(on_centre,on_radius)
40
41 
    # Define number of OFF regions for the reflected regions approach
42 
    num_off=5
43
44 
    # Define the optimizer to be used
45 
    opt = GOptimizerLM()
46
47 
    # Optionally set a spectral shape to be fitted
48 
    flux=...
49 
    index=...
50 
    spec=GModelSpectralPlaw(flux, index)
51
52 
    # Set background maker
53 
    ######################
54
55 
    # Create instance of object
56 
    maker=GOnOffMakerReflected(obs,ebds,on_reg,num_off)
57
58 
    # Compute ON-OFF data
59 
    # (loop on observation runs and energy bins to sum counts)
60 
    maker.compute_data()
61
62 
    # We can inspect the resulting array of GOnOffBin elements
63 
    for bin in maker.data():
64 
            print 'Number of ON counts: %d\n' % (bin.oncts())
65 
            print 'Number of OFF counts: %d\n' % (bin.offcts())
66 
            print 'Alpha parameter: %f\n' % (bin.alpha())
67
68 
    # Set ON-OFF fitter for two types of analysis: flux points determination or spectral shape fitting
69 
    ##################################################################################################
70
71 
    # Create instance of spectral analysis object
72 
    fitter_points= GOnOffFitterSpectrum(obs,maker)
73 
    fitter_shape= GOnOffFitterSpectrum(obs,maker,spec)
74
75 
    # Define if analysis is done on summed ON and OFF data or if run information is retained
76 
    # (this sets an internal flag)
77 
    fitter_points.sum_data()
78 
    fitter_shape.sum_data()
79
80 
    # Preparatory steps
81 
    # (compute number of expected counts for an initial spectrum, and determine true energies)
82 
    # (note that for flux points determination, an initial spectral shape is created internally)
83 
    # (these may be grouped and hidden in a prepare() method, or automatically done on instantiation ?)
84 
    fitter_points.convolve_spectrum()
85 
    fitter_points.compute_energy()
86 
    fitter_shape.convolve_spectrum()
87 
    fitter_shape.compute_energy()
88
89 
    # Perform analysis
90 
    ##################
91
92 
    # Compute/fit the best source parameters
93 
    # (implies internal calls to two different optimization functions for both analyses)
94 
    fitter_points.optimize(opt)
95 
    fitter_shape.optimize(opt)
96
97 
    # Compute significances and print them for all energy bins
98 
    fitter_points.compute_significance()
99 
    for nrj,sig in zip(fitter_points.true_energy(),fitter_points.significance()):
100 
            print "Significance %e at true energy %e\n" % (nrj,sig)
101
102 
    # Print best-fit spectral parameters
103 
    print fitter_shape.spectral()['Prefactor'].value(),' +/- ',fitter_shape.spectral()['Prefactor'].error()
104 
    print fitter_shape.spectral()['Index'].value(),' +/- ',fitter_shape.spectral()['Index'].error()
105
106 
    # Then the best-fit spectral models can be written to XML elements and then to an XML model file...
107
108
109
110
111
112