Justine Devin / gammalib

Commit 1fad59dd360720aa7c721c46a39a1317ff7d41c2

Authored by Pierrick Martin
1 parent 8b023bd1

Solved the problem in fitting. ON-OFF analysis works now.

Inline Side-by-side
Showing 1 changed file with 42 additions and 38 deletions
inst/cta/src/GCTAOnOffObservation.cpp
  View file @1fad59d
... ... @@ -691,6 +691,7 @@ double GCTAOnOffObservation::model_on(const GModels& models,
691 691 const GEnergy emin=m_on_spec.ebounds().emin(ibin);
692 692 const GEnergy emax=m_on_spec.ebounds().emax(ibin);
693 693 const GEnergy emean=m_on_spec.ebounds().elogmean(ibin);
  694 + const double ewidth=m_on_spec.ebounds().ewidth(ibin).MeV();
694 695  
695 696 // Loop over models
696 697 for (int j = 0; j < models.size(); ++j) {
... ... @@ -734,7 +735,7 @@ double GCTAOnOffObservation::model_on(const GModels&amp; models,
734 735 for (int k = 0; k < speskyptr->size(); ++k) {
735 736 GModelPar sppar=(*speskyptr)[k];
736 737 if (sppar.is_free() && i_par < npars) {
737   - (*mod_grad)[i_par]=sppar.gradient();
  738 + (*mod_grad)[i_par]=sppar.factor_gradient()*m_arf[ibin]*m_ontime*ewidth;
738 739 i_par++;
739 740 }
740 741  
... ... @@ -865,7 +866,7 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
865 866 for (int k = 0; k < spebgdptr->size(); ++k) {
866 867 GModelPar sppar=(*spebgdptr)[k];
867 868 if (sppar.is_free() && i_par < npars) {
868   - (*mod_grad)[i_par]=sppar.gradient()*m_offtime*totsolidangle*ewidth;
  869 + (*mod_grad)[i_par]=sppar.factor_gradient()*m_offtime*totsolidangle*ewidth;
869 870 i_par++;
870 871 }
871 872  
... ... @@ -915,7 +916,7 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
915 916 for (int k = 0; k < spebgdptr->size(); ++k) {
916 917 GModelPar sppar=(*spebgdptr)[k];
917 918 if (sppar.is_free() && i_par < npars) {
918   - (*mod_grad)[i_par]=sppar.gradient()*m_offtime*totsolidangle*ewidth;
  919 + (*mod_grad)[i_par]=sppar.factor_gradient()*m_offtime*totsolidangle*ewidth;
919 920 i_par++;
920 921 }
921 922  
... ... @@ -1010,8 +1011,7 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1010 1011 GVector sky_grad(npars);
1011 1012 GVector bgd_grad(npars);
1012 1013 // Working arrays
1013   - double* values = new double[npars];
1014   - int* inx = new int[npars];
  1014 + GVector colvar(npars);
1015 1015  
1016 1016 // Create time object (empty, just needed in some calls)
1017 1017 GTime time;
... ... @@ -1044,9 +1044,16 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1044 1044 nbgd=model_off(models,pars,i,&bgd_grad);
1045 1045 ngam=model_on(models,pars,i,&sky_grad);
1046 1046  
1047   - // For debug (to be removed later)
1048   - std::cout << "Obs " << m_name << " Bin " << i << " Non=" << non << " Noff=" << noff << " Ngam=" << ngam << " Nbgd=" << nbgd << std::endl;
1049   -
  1047 + // Debug: print gradient (to be removed later)
  1048 + /*
  1049 + std::cout << "Obs " << m_name << " Bin " << i << " Non=" << non << " Noff=" << noff << " Ngam=" << ngam << " Nbgd=" << nbgd << std::endl;
  1050 + for (int j = 0; j < npars; ++j) {
  1051 + if (sky_grad[j] != 0.0 || bgd_grad[j] != 0.0) {
  1052 + std::cout << " Gradient " << j << " " << sky_grad[j] << " " << bgd_grad[j] << " " << (pars[j])->name() << std::endl;
  1053 + }
  1054 + }
  1055 + */
  1056 +
1050 1057 // Skip bin if model is too small (avoids -Inf or NaN gradients)
1051 1058 nonpred= ngam+m_alpha*nbgd;
1052 1059 if ((nbgd <= minmod) || (nonpred <= minmod)) {
... ... @@ -1067,17 +1074,6 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1067 1074 sum_data += non;
1068 1075 sum_model += nonpred;
1069 1076 #endif
1070   -
1071   - // Create index array of non-zero derivatives and initialise working array
1072   - // (just needed in call to update hessian matrix)
1073   - int ndev = 0;
1074   - for (int j = 0; j < npars; ++j) {
1075   - if ((sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) ||
1076   - (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j]))) {
1077   - inx[ndev] = j;
1078   - ndev++;
1079   - }
1080   - }
1081 1077  
1082 1078 // Fill derivatives
1083 1079 double fa=non/nonpred;
... ... @@ -1100,21 +1096,22 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1100 1096  
1101 1097 // If spectral model for sky component is non-zero and non-infinite
1102 1098 if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
1103   - values[k]=sky_grad[j]*sky_grad[k]*fb;
  1099 + colvar[k]=sky_grad[j]*sky_grad[k]*fb;
1104 1100  
1105 1101 // If spectral model for sky component is non-zero and non-infinite
1106 1102 } else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
1107   - values[k]=sky_grad[j]*bgd_grad[k]*fc;
  1103 + colvar[k]=sky_grad[j]*bgd_grad[k]*fc;
1108 1104  
1109 1105 // ...else neither sky nor background
1110 1106 } else {
1111   - values[k]=0.0;
  1107 + colvar[k]=0.0;
1112 1108 }
1113 1109  
1114   - // Update matrix
1115   - curvature->add_to_column(j, values, inx, ndev);
1116   -
1117 1110 }
  1111 +
  1112 + // Update matrix
  1113 + curvature->add_to_column(j, colvar);
  1114 +
1118 1115  
1119 1116 // If spectral model for sky component is non-zero and non-infinite
1120 1117 } else if (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j])) {
... ... @@ -1127,28 +1124,39 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1127 1124  
1128 1125 // If spectral model for sky component is non-zero and non-infinite
1129 1126 if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
1130   - values[k]=bgd_grad[j]*bgd_grad[k]*fc;
  1127 + colvar[k]=bgd_grad[j]*bgd_grad[k]*fc;
1131 1128  
1132 1129 // If spectral model for sky component is non-zero and non-infinite
1133 1130 } else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
1134   - values[k]=bgd_grad[j]*bgd_grad[k]*fd;
  1131 + colvar[k]=bgd_grad[j]*bgd_grad[k]*fd;
1135 1132  
1136 1133 // ...else neither sky nor background
1137 1134 } else {
1138   - values[k]=0.0;
  1135 + colvar[k]=0.0;
1139 1136 }
1140   -
1141   - // Update matrix
1142   - curvature->add_to_column(j, values, inx, ndev);
1143   -
  1137 +
1144 1138 }
  1139 +
  1140 + // Update matrix
  1141 + curvature->add_to_column(j, colvar);
1145 1142  
1146   - }
  1143 + }
1147 1144  
1148 1145 } // Looped over all parameters for derivatives computation
1149 1146  
1150 1147 } // Looped over energy bins
1151 1148  
  1149 + // Debug: print curvature (to be removed later)
  1150 + /*
  1151 + for (int j = 0; j < npars; ++j) {
  1152 + for (int k = 0; k < npars; ++k) {
  1153 + if ((*curvature)(j,k) != 0.0) {
  1154 + std::cout << " Curvature " << j << " " << k << " " << (*curvature)(j,k) << std::endl;
  1155 + }
  1156 + }
  1157 + }
  1158 + */
  1159 +
1152 1160 // ... else parameter array is not of type GModel so raise exception
1153 1161 } else {
1154 1162 std::string msg ="No model parameter for the computation of the "
... ... @@ -1156,11 +1164,7 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1156 1164 " (ID "+this->id()+").\n";
1157 1165 throw GException::invalid_value(G_POISSON_ONOFF,msg);
1158 1166 }
1159   -
1160   - // Free working array
1161   - if (values != NULL) delete [] values;
1162   - if (inx != NULL) delete [] inx;
1163   -
  1167 +
1164 1168 // Dump statistics
1165 1169 #if G_OPT_DEBUG
1166 1170 std::cout << "Number of bins: " << n_bins << std::endl;
... ...