search

VERITAS-Observatory / EventDisplay_v4

A reconstruction and analysis pipeline for VERITAS.
BSD 3-Clause "New" or "Revised" License
4 stars 1 forks source link

VLikelihoodFitter in v487 5% disagreement #164

Closed steob92 closed 12 months ago

steob92 commented 2 years ago

In v487 VLikelihoodFitter is showing a consistent ~5% disagreement with the standard chi^2 fitter (VEnergySpectrum).

crabComparison

Standard Fit:

Fitfunction: power law
 FCN=10.7727 FROM MINOS     STATUS=SUCCESSFUL     28 CALLS         169 TOTAL
                     EDM=9.38078e-11    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           3.34774e-11   1.97291e-12  -3.08387e-14  -4.57015e+07
   2  p1          -2.64438e+00   7.74968e-02   7.74968e-02  -2.90764e-03
VSpectralFitter::fit Getting Confidence Interval (cl = 0.68)
Fit Results

Results for power law fit: 
--------------------------
dN/dE = I x (E/1 TeV)^{-Gamma}
I = 3.35e-11 +- 1.97e-12 cm^-2s^-1TeV^-1
Gamma = -2.64 +- 0.08
Chi2 10.77, N = 4 (Chi2/N=2.69)

Likelihood Fit

 **********
 **    8 **MINOS       1e+04           2
 **********
 FCN=-2933.61 FROM MINOS     STATUS=SUCCESSFUL     25 CALLS         129 TOTAL
                     EDM=1.7422e-08    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  Norm         3.20604e-11   1.85889e-12  -1.83265e-12   1.88842e-12
   2  Index       -2.52458e+00   5.38055e-02  -5.43379e-02   5.33839e-02
                               ERR DEF= 0.5
Variabile 1 Error Status: 1, E_Low = -5.43e-02, E_Up = 5.34e-02
Binned Likelihood Fit: 
Parameter    Best Fit    ErrorL      ErrorU      IsMin
[0]      Norm       3.21e-11        -1.83e-12       1.89e-12        1
[1]      Index      -2.52e+00       -5.43e-02       5.34e-02        1
E_Norm: 1.00e+00

Printing Covariance Matrix
3.46e-24    7.23e-14    
7.23e-14    2.90e-03    

Calculating Total Chi^2
Total Live time (s): 4.05e+03
Number of Bins: 10
Number of Fitting Parameters: 2
Number of Degrees of Freedom: 8
Chi^2 :2.12e+01
Chi^2/NDF : 2.66e+00
log(L) : 2.93e+03
log(L_0) : 2.94e+03
steob92 commented 2 years ago

Here is the same plot from ED v486. The results are consistent with no trend. Looks like this is isolated to v487.

The runs used in this analysis are

64080
64081
64082
64083

crabComparison_v486 .

steob92 commented 2 years ago

Standard Fit Results (v486):

Fitfunction: power law
 FCN=10.7727 FROM MINOS     STATUS=SUCCESSFUL     28 CALLS         169 TOTAL
                     EDM=9.38078e-11    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           3.34774e-11   1.97291e-12  -3.08387e-14  -4.57015e+07
   2  p1          -2.64438e+00   7.74968e-02   7.74968e-02  -2.90764e-03
VSpectralFitter::fit Getting Confidence Interval (cl = 0.68)
Fit Results

Results for power law fit: 
--------------------------
dN/dE = I x (E/1 TeV)^{-Gamma}
I = 3.35e-11 +- 1.97e-12 cm^-2s^-1TeV^-1
Gamma = -2.64 +- 0.08
Chi2 10.77, N = 4 (Chi2/N=2.69)

Likelihood Fit Results (v486):

Variabile 0 Error Status: 1, E_Low = -1.94e-12, E_Up = 2.00e-12
 **********
 **    8 **MINOS       1e+04           2
 **********
 FCN=-2933.71 FROM MINOS     STATUS=SUCCESSFUL     25 CALLS         122 TOTAL
                     EDM=6.99485e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  Norm         3.37660e-11   1.96874e-12  -1.93993e-12   2.00119e-12
   2  Index       -2.54355e+00   5.43363e-02  -5.48862e-02   5.38996e-02
                               ERR DEF= 0.5
Variabile 1 Error Status: 1, E_Low = -5.49e-02, E_Up = 5.39e-02
Binned Likelihood Fit: 
Parameter    Best Fit    ErrorL      ErrorU      IsMin
[0]      Norm       3.38e-11        -1.94e-12       2.00e-12        1
[1]      Index      -2.54e+00       -5.49e-02       5.39e-02        1
E_Norm: 1.00e+00

Printing Covariance Matrix
3.88e-24    7.77e-14    
7.77e-14    2.95e-03    

Calculating Total Chi^2
Total Live time (s): 4.05e+03
Number of Bins: 10
Number of Fitting Parameters: 2
Number of Degrees of Freedom: 8
Chi^2 :2.10e+01
Chi^2/NDF : 2.63e+00
log(L) : 2.93e+03
log(L_0) : 2.94e+03
GernotMaier commented 2 years ago

OK. That will be a nightmare to debug.

IRFs should be exactly the same between v487 and v486 (with the exception of the axis bug in e_sys_rel). There is no new IRF calculation for v487, it is only filled differently at the combine effective area stage. This requires also a different reading of the values in VEffectiveArea.

GernotMaier commented 2 years ago

Just checked the gammapy analysis between v486 and v487, relying on the same IRFs as the binned likelihood fitter in anasum: results are identical (<0.05% difference). This means it is an issue with reading the values.

GernotMaier commented 2 years ago

@steob92 please test again with the last commit c74a45b

It was indeed an issue with reading the values: with two effective areas trees (one for values as function of reconstructed energies and the second one as function of true energies), one needs to take care to read out the correct values.

effective area tree with 20x less entries and the relation between the tree entry index is:

index_eff_true = index_eff_rec / n_index_bins + index_eff_rec % n_az_bins

VEffectiveAreaCalculator is incredibly complicated.

steob92 commented 2 years ago

I'll take a look. Is the previous effective area file you game me still valid?

GernotMaier commented 2 years ago

Yes, that file should work.

I was simply getting the wrong entry from the tree for the IRFs as function of true energy.

steob92 commented 2 years ago

The results look good.

There is a very small (~2%) discrepancy in the lowest energy bins. This isn't new, but its a little more evident than the v486 analysis. I believe it is due to how threshold are handled in standard analysis and likelihood analysis. The response matrix essentially includes sub-threshold events when integrating to get the predicted above threshold events. I don't know if this is worth pursuing more?

Standard Fit results (v487):

Fitfunction: power law
 FCN=12.7193 FROM MINOS     STATUS=SUCCESSFUL     28 CALLS         157 TOTAL
                     EDM=2.46705e-08    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           3.32552e-11   1.87543e-12  -3.26915e-14  -5.18097e+07
   2  p1          -2.65445e+00   7.13930e-02   7.13930e-02  -1.10453e-04
VSpectralFitter::fit Getting Confidence Interval (cl = 0.7)
IGNORING FITTER
Fit Results

Results for power law fit: 
--------------------------
dN/dE = I x (E/1.0 TeV)^{-Gamma}
I = 3.33e-11 +- 1.88e-12 cm^-2s^-1TeV^-1
Gamma = -2.65 +- 0.07
Chi2 12.72, N = 6 (Chi2/N=2.12)

Likelihood Fit result:

Variabile 0 Error Status: 1, E_Low = -1.97e-12, E_Up = 2.03e-12
 **********
 **    8 **MINOS       1e+04           2
 **********
 FCN=-2933.71 FROM MINOS     STATUS=SUCCESSFUL     25 CALLS         122 TOTAL
                     EDM=4.76278e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  Norm         3.42573e-11   1.99727e-12  -1.97051e-12   2.02750e-12
   2  Index       -2.54677e+00   5.40599e-02  -5.46231e-02   5.36125e-02
                               ERR DEF= 0.5
Variabile 1 Error Status: 1, E_Low = -5.46e-02, E_Up = 5.36e-02
Binned Likelihood Fit: 
Parameter    Best Fit    ErrorL      ErrorU      IsMin
[0]      Norm       3.43e-11        -1.97e-12       2.03e-12        1
[1]      Index      -2.55e+00       -5.46e-02       5.36e-02        1
E_Norm: 1.00e+00

Printing Covariance Matrix
3.99e-24    7.85e-14    
7.85e-14    2.92e-03    

Calculating Total Chi^2
Total Live time (s): 4.05e+03
Number of Bins: 10
Number of Fitting Parameters: 2
Number of Degrees of Freedom: 8
Chi^2 :2.11e+01
Chi^2/NDF : 2.63e+00
log(L) : 2.93e+03
log(L_0) : 2.94e+03

crabComparison

ratio

GernotMaier commented 2 years ago

But this is a difference between standard and likelihood analysis, and not a difference between the two different ways of writing IRFs?

steob92 commented 2 years ago

Correct. This is unrelated to how the IRFs are generated.

GernotMaier commented 2 years ago

Good - I am closing the issue, as I think we have solved the initial problem.

steob92 commented 2 years ago

Hi Gernot,

Was this bug fix implemented in v487b? I've been looking into the previously mention 2% discrepancy around energy threshold and it appears that I'm seeing a consistent 5% discrepancy above ~700 GeV. We previously saw this in a development branch for v487 (this issue).

Before I go into further debugging, was this fix implemented into v487b?

All the best, Ste

CrabCompare ratio

GernotMaier commented 2 years ago

Just looked again through the code and the fix is applied (was commit c74a45bda1f97c90e04e68d127fafa8493430e51).

It is already fixed in v487b (see https://github.com/VERITAS-Observatory/EventDisplay_v4/blob/87bd52c5b769d6f4c42c0194c7ea80fb6211e439/src/VEffectiveAreaCalculator.cpp#L1369).

Not sure what to do now - are you absolutely sure that it re-appeared or was it never gone?

steob92 commented 2 years ago

I'm not entirely sure I'm not just doing something stupid...

I have some local changes that I was testing. I'm also using IRFs downloaded by someone else. Let me try and consolidate with a base v487b install and IRFs I can confirm are the correct versions. I'll take a look over the weekend.

GernotMaier commented 2 years ago

Let's reopen this issue and address it. Looks like Sonal is seeing something similar when looking into the gammapy results.

steob92 commented 2 years ago

Rerunning knowing I'm using the correct IRFs and using a fresh v487c install I still see the same results for the previously posted runlist, which was:

RUN 69976 (+0.0N,-0.5W)  at 41,  87 deg El., Az, 30.05 min, Non:  104, Noff: 10.67 (  64, norm 0.167),  14.2 sigma, Rates:   3.106 +/-   0.342 gamma/min (background:   0.355 events/min)
RUN 74506 (+0.0N,-0.5W)  at 80, 172 deg El., Az, 20.02 min, Non:  166, Noff:  6.17 (  37, norm 0.167),  21.6 sigma, Rates:   7.985 +/-   0.646 gamma/min (background:   0.308 events/min)
RUN 74860 (-0.5N,-0.0W)  at 65,-109 deg El., Az,  8.42 min, Non:   62, Noff:  2.50 (  15, norm 0.167),  13.0 sigma, Rates:   7.069 +/-   0.939 gamma/min (background:   0.297 events/min)
RUN 80926 (-0.0N,+0.5W)  at 54, -96 deg El., Az, 30.02 min, Non:  179, Noff:  5.33 (  32, norm 0.167),  23.0 sigma, Rates:   5.786 +/-   0.447 gamma/min (background:   0.178 events/min)
RUN 84175 (+0.5N,+0.0W)  at 47,  90 deg El., Az, 30.02 min, Non:  110, Noff:  5.50 (  33, norm 0.167),  16.8 sigma, Rates:   3.481 +/-   0.351 gamma/min (background:   0.183 events/min)

Using the Crab_V6_v400_release runlist:

64080
64081
64082
64083

I get the attached results. The comparison is a little better, but there still is a ~2.5% discrepancy.

CrabCompare ratio

I've also attached the MC effective areas for run 64081 (RUN 64081 (-0.5N,-0.0W) at 79, 170 deg El., Az, 20.03 min) for v483, v485 and v487 and the ratio of v4XX/v487.

Eff_64081_ratio Eff_64081

GernotMaier commented 2 years ago

@steob92 - I think I found the issue:

First, I can replicate the problem (there is an option when combining effective area trees (DL3test) which writes the IRFs vs true energy both in the fEffArea and fEffAreaH2F truee and allows testing. Here an example of discrepancies of effective area values:

Screenshot 2022-05-31 at 20 50 17

I bit of digging found hopefully the issue and I've applied a bug fix in the branch v487d-dev-v0.1. This is a single line which changed, see here

Could you check again and test the problem with v487d-dev-v0.1?

steob92 commented 2 years ago

Unfortunatly this doesn't appear to solve the issue. The discrepancy increases to ~5% using the Crab_V6_v400_release runlist:

64080
64081
64082
64083

CrabCompare ratio

GernotMaier commented 12 months ago

I am closing this, as we do not expect to make progress on this topic