Changeset 8911 for trunk/MagicSoft


Ignore:
Timestamp:
06/02/08 10:36:33 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc

    r8359 r8911  
    1818!   Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2004
     20!   Copyright: MAGIC Software Development, 2000-2008
    2121!
    2222!
    … …  
    432432    {
    433433    case kEPoisson4:
    434       fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
     434      fSinglePheFit = new TF1("SinglePheFit",&PoissonKto4,rmin,rmax,6);
    435435      rmin += 6.5;
    436436      break;
    437437    case kEPoisson5:
    438       fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
     438      fSinglePheFit = new TF1("SinglePheFit",&PoissonKto5,rmin,rmax,6);
    439439      rmin = 0.;
    440440      break;
    441441    case kEPoisson6:
    442       fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
     442      fSinglePheFit = new TF1("SinglePheFit",&PoissonKto6,rmin,rmax,6);
    443443      break;
    444444    case kEPolya:
    445       fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
     445      fSinglePheFit = new TF1("SinglePheFit",&Polya,rmin,rmax,8);
    446446      break;
    447447    case kEMichele:
    448       fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9);
     448      fSinglePheFit = new TF1("SinglePheFit",&FitFuncMichele,rmin,rmax,9);
    449449      break;
    450450    default:
    … …  
    964964}
    965965
    966 
    967 
    968 
    969 
    970 
    971 
    972 
    973 
    974 
    975 
     966Double_t MHCalibrationChargeBlindPix::FitFuncMichele(Double_t *x, Double_t *par)
     967{
     968    Double_t lambda1cat = par[0];
     969    Double_t lambda1dyn = par[1];
     970    Double_t mu0        = par[2];
     971    Double_t mu1cat     = par[3];
     972    Double_t mu1dyn     = par[4];
     973    Double_t sigma0     = par[5];
     974    Double_t sigma1cat  = par[6];
     975    Double_t sigma1dyn  = par[7];
     976
     977    Double_t sumcat = 0.;
     978    Double_t sumdyn = 0.;
     979    Double_t arg    = 0.;
     980
     981    if (lambda1cat < lambda1dyn)
     982        return FLT_MAX;
     983
     984    if (mu1cat    < mu0)
     985        return FLT_MAX;
     986
     987    if (mu1dyn    < mu0)
     988        return FLT_MAX;
     989
     990    if (mu1cat < mu1dyn)
     991        return FLT_MAX;
     992
     993    if (sigma0 < 0.0001)
     994        return FLT_MAX;
     995
     996    if (sigma1cat < sigma0)
     997        return FLT_MAX;
     998
     999    if (sigma1dyn < sigma0)
     1000        return FLT_MAX;
     1001
     1002    Double_t mu2cat = (2.*mu1cat)-mu0;
     1003    Double_t mu2dyn = (2.*mu1dyn)-mu0;
     1004    Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
     1005    Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
     1006
     1007    Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0));
     1008    Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0));
     1009    Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
     1010    Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
     1011
     1012    Double_t lambda2cat = lambda1cat*lambda1cat;
     1013    Double_t lambda2dyn = lambda1dyn*lambda1dyn;
     1014    Double_t lambda3cat = lambda2cat*lambda1cat;
     1015    Double_t lambda3dyn = lambda2dyn*lambda1dyn;
     1016
     1017    // k=0:
     1018    arg = (x[0] - mu0)/sigma0;
     1019    sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
     1020    sumdyn = sumcat;
     1021
     1022    // k=1cat:
     1023    arg = (x[0] - mu1cat)/sigma1cat;
     1024    sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
     1025    // k=1dyn:
     1026    arg = (x[0] - mu1dyn)/sigma1dyn;
     1027    sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
     1028
     1029    // k=2cat:
     1030    arg = (x[0] - mu2cat)/sigma2cat;
     1031    sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
     1032    // k=2dyn:
     1033    arg = (x[0] - mu2dyn)/sigma2dyn;
     1034    sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
     1035
     1036
     1037    // k=3cat:
     1038    arg = (x[0] - mu3cat)/sigma3cat;
     1039    sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
     1040    // k=3dyn:
     1041    arg = (x[0] - mu3dyn)/sigma3dyn;
     1042    sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn;
     1043
     1044    sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
     1045    sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
     1046
     1047    return par[8]*(sumcat+sumdyn)/2.;
     1048}
     1049
     1050Double_t MHCalibrationChargeBlindPix::PoissonKto4(Double_t *x, Double_t *par)
     1051{
     1052    Double_t lambda = par[0];
     1053
     1054    Double_t sum = 0.;
     1055    Double_t arg = 0.;
     1056
     1057    Double_t mu0 = par[1];
     1058    Double_t mu1 = par[2];
     1059
     1060    if (mu1 < mu0)
     1061        return FLT_MAX;
     1062
     1063    Double_t sigma0 = par[3];
     1064    Double_t sigma1 = par[4];
     1065
     1066    if (sigma0 < 0.0001)
     1067        return FLT_MAX;
     1068
     1069    if (sigma1 < sigma0)
     1070        return FLT_MAX;
     1071
     1072    Double_t mu2 = (2.*mu1)-mu0;
     1073    Double_t mu3 = (3.*mu1)-(2.*mu0);
     1074    Double_t mu4 = (4.*mu1)-(3.*mu0);
     1075
     1076    Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
     1077    Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
     1078    Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
     1079
     1080    Double_t lambda2 = lambda*lambda;
     1081    Double_t lambda3 = lambda2*lambda;
     1082    Double_t lambda4 = lambda3*lambda;
     1083
     1084    // k=0:
     1085    arg = (x[0] - mu0)/sigma0;
     1086    sum = TMath::Exp(-0.5*arg*arg)/sigma0;
     1087
     1088    // k=1:
     1089    arg = (x[0] - mu1)/sigma1;
     1090    sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
     1091
     1092    // k=2:
     1093    arg = (x[0] - mu2)/sigma2;
     1094    sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
     1095
     1096    // k=3:
     1097    arg = (x[0] - mu3)/sigma3;
     1098    sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
     1099
     1100    // k=4:
     1101    arg = (x[0] - mu4)/sigma4;
     1102    sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
     1103
     1104    return TMath::Exp(-1.*lambda)*par[5]*sum;
     1105}
     1106
     1107Double_t MHCalibrationChargeBlindPix::PoissonKto5(Double_t *x, Double_t *par)
     1108{
     1109    Double_t lambda = par[0];
     1110
     1111    Double_t sum = 0.;
     1112    Double_t arg = 0.;
     1113
     1114    Double_t mu0 = par[1];
     1115    Double_t mu1 = par[2];
     1116
     1117    if (mu1 < mu0)
     1118        return FLT_MAX;
     1119
     1120    Double_t sigma0 = par[3];
     1121    Double_t sigma1 = par[4];
     1122
     1123    if (sigma0 < 0.0001)
     1124        return FLT_MAX;
     1125
     1126    if (sigma1 < sigma0)
     1127        return FLT_MAX;
     1128
     1129
     1130    Double_t mu2 = (2.*mu1)-mu0;
     1131    Double_t mu3 = (3.*mu1)-(2.*mu0);
     1132    Double_t mu4 = (4.*mu1)-(3.*mu0);
     1133    Double_t mu5 = (5.*mu1)-(4.*mu0);
     1134
     1135    Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
     1136    Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
     1137    Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
     1138    Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
     1139
     1140    Double_t lambda2 = lambda*lambda;
     1141    Double_t lambda3 = lambda2*lambda;
     1142    Double_t lambda4 = lambda3*lambda;
     1143    Double_t lambda5 = lambda4*lambda;
     1144
     1145    // k=0:
     1146    arg = (x[0] - mu0)/sigma0;
     1147    sum = TMath::Exp(-0.5*arg*arg)/sigma0;
     1148
     1149    // k=1:
     1150    arg = (x[0] - mu1)/sigma1;
     1151    sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
     1152
     1153    // k=2:
     1154    arg = (x[0] - mu2)/sigma2;
     1155    sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
     1156
     1157    // k=3:
     1158    arg = (x[0] - mu3)/sigma3;
     1159    sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
     1160
     1161    // k=4:
     1162    arg = (x[0] - mu4)/sigma4;
     1163    sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
     1164
     1165    // k=5:
     1166    arg = (x[0] - mu5)/sigma5;
     1167    sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
     1168
     1169    return TMath::Exp(-1.*lambda)*par[5]*sum;
     1170}
     1171
     1172Double_t MHCalibrationChargeBlindPix::PoissonKto6(Double_t *x, Double_t *par)
     1173{
     1174    Double_t lambda = par[0];
     1175
     1176    Double_t sum = 0.;
     1177    Double_t arg = 0.;
     1178
     1179    Double_t mu0 = par[1];
     1180    Double_t mu1 = par[2];
     1181
     1182    if (mu1 < mu0)
     1183        return FLT_MAX;
     1184
     1185    Double_t sigma0 = par[3];
     1186    Double_t sigma1 = par[4];
     1187
     1188    if (sigma0 < 0.0001)
     1189        return FLT_MAX;
     1190
     1191    if (sigma1 < sigma0)
     1192        return FLT_MAX;
     1193
     1194
     1195    Double_t mu2 = (2.*mu1)-mu0;
     1196    Double_t mu3 = (3.*mu1)-(2.*mu0);
     1197    Double_t mu4 = (4.*mu1)-(3.*mu0);
     1198    Double_t mu5 = (5.*mu1)-(4.*mu0);
     1199    Double_t mu6 = (6.*mu1)-(5.*mu0);
     1200
     1201    Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
     1202    Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
     1203    Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
     1204    Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
     1205    Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
     1206
     1207    Double_t lambda2 = lambda*lambda;
     1208    Double_t lambda3 = lambda2*lambda;
     1209    Double_t lambda4 = lambda3*lambda;
     1210    Double_t lambda5 = lambda4*lambda;
     1211    Double_t lambda6 = lambda5*lambda;
     1212
     1213    // k=0:
     1214    arg = (x[0] - mu0)/sigma0;
     1215    sum = TMath::Exp(-0.5*arg*arg)/sigma0;
     1216
     1217    // k=1:
     1218    arg = (x[0] - mu1)/sigma1;
     1219    sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
     1220
     1221    // k=2:
     1222    arg = (x[0] - mu2)/sigma2;
     1223    sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
     1224
     1225    // k=3:
     1226    arg = (x[0] - mu3)/sigma3;
     1227    sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
     1228
     1229    // k=4:
     1230    arg = (x[0] - mu4)/sigma4;
     1231    sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
     1232
     1233    // k=5:
     1234    arg = (x[0] - mu5)/sigma5;
     1235    sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
     1236
     1237    // k=6:
     1238    arg = (x[0] - mu6)/sigma6;
     1239    sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
     1240
     1241    return TMath::Exp(-1.*lambda)*par[5]*sum;
     1242}
     1243
     1244Double_t MHCalibrationChargeBlindPix::Polya(Double_t *x, Double_t *par)
     1245{
     1246    const Double_t QEcat = 0.247;            // mean quantum efficiency
     1247    const Double_t sqrt2 = 1.4142135623731;
     1248    const Double_t sqrt3 = 1.7320508075689;
     1249    const Double_t sqrt4 = 2.;
     1250
     1251    const Double_t lambda = par[0];           // mean number of photons
     1252
     1253    const Double_t excessPoisson = par[1];    // non-Poissonic noise contribution
     1254    const Double_t delta1 = par[2];           // amplification first dynode
     1255    const Double_t delta2 = par[3];           // amplification subsequent dynodes
     1256
     1257    const Double_t electronicAmpl = par[4];   // electronic amplification and conversion to FADC charges
     1258
     1259    const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2;  // total PMT gain
     1260    const Double_t A = 1. + excessPoisson - QEcat
     1261        + 1./delta1
     1262        + 1./delta1/delta2
     1263        + 1./delta1/delta2/delta2;                                  // variance contributions from PMT and QE
     1264
     1265    const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl;        // Total gain and conversion
     1266
     1267    const Double_t mu0 = par[7];                                      // pedestal
     1268    const Double_t mu1 = totAmpl;                                 // single phe position
     1269    const Double_t mu2 = 2*totAmpl;                               // double phe position
     1270    const Double_t mu3 = 3*totAmpl;                               // triple phe position
     1271    const Double_t mu4 = 4*totAmpl;                               // quadruple phe position
     1272
     1273    const Double_t sigma0 = par[5];
     1274    const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
     1275    const Double_t sigma2 = sqrt2*sigma1;
     1276    const Double_t sigma3 = sqrt3*sigma1;
     1277    const Double_t sigma4 = sqrt4*sigma1;
     1278
     1279    const Double_t lambda2 = lambda*lambda;
     1280    const Double_t lambda3 = lambda2*lambda;
     1281    const Double_t lambda4 = lambda3*lambda;
     1282
     1283    //-- calculate the area----
     1284    Double_t arg = (x[0] - mu0)/sigma0;
     1285    Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
     1286
     1287    // k=1:
     1288    arg = (x[0] - mu1)/sigma1;
     1289    sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
     1290
     1291    // k=2:
     1292    arg = (x[0] - mu2)/sigma2;
     1293    sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
     1294
     1295    // k=3:
     1296    arg = (x[0] - mu3)/sigma3;
     1297    sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
     1298
     1299    // k=4:
     1300    arg = (x[0] - mu4)/sigma4;
     1301    sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
     1302
     1303    return TMath::Exp(-1.*lambda)*par[6]*sum;
     1304}

  • trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h

    r8359 r8911  
    55#include "MHCalibrationPix.h"
    66#endif
     7
     8#include <TMath.h>
    79
    810#ifndef ROOT_TF1
    … …  
    5961  Byte_t    fFlags;                          // Bit-field for the flags
    6062  enum { kSinglePheFitOK, kPedestalFitOK };  // Possible bits to be set
     63
     64  static Double_t FitFuncMichele(Double_t *x, Double_t *par);
     65  static Double_t PoissonKto4(Double_t *x, Double_t *par);
     66  static Double_t PoissonKto5(Double_t *x, Double_t *par);
     67  static Double_t PoissonKto6(Double_t *x, Double_t *par);
     68  static Double_t Polya(Double_t *x, Double_t *par);
    6169
    6270public:
    … …  
    137145                           const Double_t sigma0, const Double_t sigma1);
    138146 
    139 private:
    140 
    141   inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par)
    142     {
    143      
    144       Double_t lambda1cat = par[0]; 
    145       Double_t lambda1dyn = par[1];
    146       Double_t mu0        = par[2];
    147       Double_t mu1cat     = par[3];
    148       Double_t mu1dyn     = par[4];
    149       Double_t sigma0     = par[5];
    150       Double_t sigma1cat  = par[6];
    151       Double_t sigma1dyn  = par[7];
    152      
    153       Double_t sumcat = 0.;
    154       Double_t sumdyn = 0.;
    155       Double_t arg    = 0.;
    156      
    157       if (lambda1cat < lambda1dyn)
    158         return FLT_MAX;
    159 
    160       if (mu1cat    < mu0)
    161         return FLT_MAX;
    162 
    163       if (mu1dyn    < mu0)
    164         return FLT_MAX;
    165 
    166       if (mu1cat < mu1dyn)
    167         return FLT_MAX;
    168 
    169       if (sigma0 < 0.0001)
    170         return FLT_MAX;
    171      
    172       if (sigma1cat < sigma0)
    173         return FLT_MAX;
    174 
    175       if (sigma1dyn < sigma0)
    176         return FLT_MAX;
    177 
    178       Double_t mu2cat = (2.*mu1cat)-mu0; 
    179       Double_t mu2dyn = (2.*mu1dyn)-mu0; 
    180       Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
    181       Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
    182      
    183       Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); 
    184       Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); 
    185       Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
    186       Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
    187      
    188       Double_t lambda2cat = lambda1cat*lambda1cat;
    189       Double_t lambda2dyn = lambda1dyn*lambda1dyn;
    190       Double_t lambda3cat = lambda2cat*lambda1cat;
    191       Double_t lambda3dyn = lambda2dyn*lambda1dyn;
    192 
    193      // k=0:
    194       arg = (x[0] - mu0)/sigma0;
    195       sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
    196       sumdyn = sumcat;
    197 
    198       // k=1cat:
    199       arg = (x[0] - mu1cat)/sigma1cat;
    200       sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
    201       // k=1dyn:
    202       arg = (x[0] - mu1dyn)/sigma1dyn;
    203       sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
    204      
    205       // k=2cat:
    206       arg = (x[0] - mu2cat)/sigma2cat;
    207       sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
    208       // k=2dyn:
    209       arg = (x[0] - mu2dyn)/sigma2dyn;
    210       sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
    211  
    212      
    213       // k=3cat:
    214       arg = (x[0] - mu3cat)/sigma3cat;
    215       sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
    216       // k=3dyn:
    217       arg = (x[0] - mu3dyn)/sigma3dyn;
    218       sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; 
    219    
    220       sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
    221       sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
    222      
    223       return par[8]*(sumcat+sumdyn)/2.;
    224 
    225     }
    226    
    227   inline static Double_t fPoissonKto4(Double_t *x, Double_t *par)
    228     {
    229 
    230       Double_t lambda = par[0]; 
    231      
    232       Double_t sum = 0.;
    233       Double_t arg = 0.;
    234      
    235       Double_t mu0 = par[1];
    236       Double_t mu1 = par[2];
    237      
    238       if (mu1 < mu0)
    239         return FLT_MAX;
    240 
    241       Double_t sigma0 = par[3];
    242       Double_t sigma1 = par[4];
    243 
    244       if (sigma0 < 0.0001)
    245         return FLT_MAX;
    246      
    247       if (sigma1 < sigma0)
    248         return FLT_MAX;
    249      
    250       Double_t mu2 = (2.*mu1)-mu0; 
    251       Double_t mu3 = (3.*mu1)-(2.*mu0);
    252       Double_t mu4 = (4.*mu1)-(3.*mu0);
    253      
    254       Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); 
    255       Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
    256       Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
    257      
    258       Double_t lambda2 = lambda*lambda;
    259       Double_t lambda3 = lambda2*lambda;
    260       Double_t lambda4 = lambda3*lambda;
    261      
    262       // k=0:
    263       arg = (x[0] - mu0)/sigma0;
    264       sum = TMath::Exp(-0.5*arg*arg)/sigma0;
    265      
    266       // k=1:
    267       arg = (x[0] - mu1)/sigma1;
    268       sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
    269      
    270       // k=2:
    271       arg = (x[0] - mu2)/sigma2;
    272       sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
    273      
    274       // k=3:
    275       arg = (x[0] - mu3)/sigma3;
    276       sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
    277      
    278       // k=4:
    279       arg = (x[0] - mu4)/sigma4;
    280       sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
    281      
    282       return TMath::Exp(-1.*lambda)*par[5]*sum;
    283      
    284     }
    285 
    286  
    287   inline static Double_t fPoissonKto5(Double_t *x, Double_t *par)
    288     {
    289      
    290       Double_t lambda = par[0]; 
    291      
    292       Double_t sum = 0.;
    293       Double_t arg = 0.;
    294      
    295       Double_t mu0 = par[1];
    296       Double_t mu1 = par[2];
    297      
    298       if (mu1 < mu0)
    299         return FLT_MAX;
    300      
    301       Double_t sigma0 = par[3];
    302       Double_t sigma1 = par[4];
    303      
    304       if (sigma0 < 0.0001)
    305         return FLT_MAX;
    306      
    307       if (sigma1 < sigma0)
    308         return FLT_MAX;
    309      
    310      
    311       Double_t mu2 = (2.*mu1)-mu0; 
    312       Double_t mu3 = (3.*mu1)-(2.*mu0);
    313       Double_t mu4 = (4.*mu1)-(3.*mu0);
    314       Double_t mu5 = (5.*mu1)-(4.*mu0);
    315      
    316       Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); 
    317       Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
    318       Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
    319       Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
    320      
    321       Double_t lambda2 = lambda*lambda;
    322       Double_t lambda3 = lambda2*lambda;
    323       Double_t lambda4 = lambda3*lambda;
    324       Double_t lambda5 = lambda4*lambda;
    325      
    326       // k=0:
    327       arg = (x[0] - mu0)/sigma0;
    328       sum = TMath::Exp(-0.5*arg*arg)/sigma0;
    329      
    330       // k=1:
    331       arg = (x[0] - mu1)/sigma1;
    332       sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
    333      
    334       // k=2:
    335       arg = (x[0] - mu2)/sigma2;
    336       sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
    337      
    338       // k=3:
    339       arg = (x[0] - mu3)/sigma3;
    340       sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
    341      
    342       // k=4:
    343       arg = (x[0] - mu4)/sigma4;
    344       sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
    345      
    346       // k=5:
    347       arg = (x[0] - mu5)/sigma5;
    348       sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
    349      
    350       return TMath::Exp(-1.*lambda)*par[5]*sum;
    351      
    352     }
    353  
    354  
    355   inline static Double_t fPoissonKto6(Double_t *x, Double_t *par)
    356     {
    357      
    358       Double_t lambda = par[0]; 
    359      
    360       Double_t sum = 0.;
    361       Double_t arg = 0.;
    362      
    363       Double_t mu0 = par[1];
    364       Double_t mu1 = par[2];
    365      
    366       if (mu1 < mu0)
    367         return FLT_MAX;
    368      
    369       Double_t sigma0 = par[3];
    370       Double_t sigma1 = par[4];
    371      
    372       if (sigma0 < 0.0001)
    373         return FLT_MAX;
    374      
    375       if (sigma1 < sigma0)
    376         return FLT_MAX;
    377      
    378      
    379       Double_t mu2 = (2.*mu1)-mu0; 
    380       Double_t mu3 = (3.*mu1)-(2.*mu0);
    381       Double_t mu4 = (4.*mu1)-(3.*mu0);
    382       Double_t mu5 = (5.*mu1)-(4.*mu0);
    383       Double_t mu6 = (6.*mu1)-(5.*mu0);
    384      
    385       Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); 
    386       Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
    387       Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
    388       Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
    389       Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
    390      
    391       Double_t lambda2 = lambda*lambda;
    392       Double_t lambda3 = lambda2*lambda;
    393       Double_t lambda4 = lambda3*lambda;
    394       Double_t lambda5 = lambda4*lambda;
    395       Double_t lambda6 = lambda5*lambda;
    396      
    397       // k=0:
    398       arg = (x[0] - mu0)/sigma0;
    399       sum = TMath::Exp(-0.5*arg*arg)/sigma0;
    400      
    401       // k=1:
    402       arg = (x[0] - mu1)/sigma1;
    403       sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
    404      
    405       // k=2:
    406       arg = (x[0] - mu2)/sigma2;
    407       sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
    408      
    409       // k=3:
    410       arg = (x[0] - mu3)/sigma3;
    411       sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
    412      
    413       // k=4:
    414       arg = (x[0] - mu4)/sigma4;
    415       sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
    416      
    417       // k=5:
    418       arg = (x[0] - mu5)/sigma5;
    419       sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
    420      
    421       // k=6:
    422       arg = (x[0] - mu6)/sigma6;
    423       sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
    424      
    425       return TMath::Exp(-1.*lambda)*par[5]*sum;
    426      
    427     }
    428 
    429   inline static Double_t fPolya(Double_t *x, Double_t *par)
    430     {
    431 
    432       const Double_t QEcat = 0.247;            // mean quantum efficiency
    433       const Double_t sqrt2 = 1.4142135623731;
    434       const Double_t sqrt3 = 1.7320508075689;
    435       const Double_t sqrt4 = 2.;
    436      
    437       const Double_t lambda = par[0];           // mean number of photons
    438      
    439       const Double_t excessPoisson = par[1];    // non-Poissonic noise contribution
    440       const Double_t delta1 = par[2];           // amplification first dynode
    441       const Double_t delta2 = par[3];           // amplification subsequent dynodes
    442      
    443       const Double_t electronicAmpl = par[4];   // electronic amplification and conversion to FADC charges
    444 
    445       const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2;  // total PMT gain
    446       const Double_t A = 1. + excessPoisson - QEcat                       
    447         + 1./delta1
    448                 + 1./delta1/delta2
    449         + 1./delta1/delta2/delta2;                                  // variance contributions from PMT and QE
    450      
    451       const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl;        // Total gain and conversion
    452      
    453       const Double_t mu0 = par[7];                                      // pedestal
    454       const Double_t mu1 = totAmpl;                                 // single phe position
    455       const Double_t mu2 = 2*totAmpl;                               // double phe position
    456       const Double_t mu3 = 3*totAmpl;                               // triple phe position
    457       const Double_t mu4 = 4*totAmpl;                               // quadruple phe position
    458      
    459       const Double_t sigma0 = par[5];
    460       const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
    461       const Double_t sigma2 = sqrt2*sigma1;
    462       const Double_t sigma3 = sqrt3*sigma1;
    463       const Double_t sigma4 = sqrt4*sigma1;
    464      
    465       const Double_t lambda2 = lambda*lambda;
    466       const Double_t lambda3 = lambda2*lambda;
    467       const Double_t lambda4 = lambda3*lambda;
    468      
    469       //-- calculate the area----
    470       Double_t arg = (x[0] - mu0)/sigma0;
    471       Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
    472      
    473      // k=1:
    474       arg = (x[0] - mu1)/sigma1;
    475       sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
    476      
    477       // k=2:
    478       arg = (x[0] - mu2)/sigma2;
    479       sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
    480      
    481       // k=3:
    482       arg = (x[0] - mu3)/sigma3;
    483       sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
    484      
    485       // k=4:
    486       arg = (x[0] - mu4)/sigma4;
    487       sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;     
    488      
    489       return TMath::Exp(-1.*lambda)*par[6]*sum;
    490     }
    491  
    492147  ClassDef(MHCalibrationChargeBlindPix, 1)  // Histogram class for Charge Blind Pixel Calibration
    493148};

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

    r8709 r8911  
    1919!   Author(s): Harald Kornmayer 1/2001
    2020!
    21 !   Copyright: MAGIC Software Development, 2000-2006
     21!   Copyright: MAGIC Software Development, 2000-2008
    2222!
    2323!
    … …  
    3535//////////////////////////////////////////////////////////////////////////////
    3636#include "MHCollectionArea.h"
     37
     38#include <TMath.h>
    3739
    3840#include <TLatex.h>
    … …  
    112114    fHistAll.Sumw2();
    113115    fHEnergy.Sumw2();
     116}
     117
     118// --------------------------------------------------------------------------
     119//
     120// Return the Area defined by fMcAreaRadius
     121//
     122Double_t MHCollectionArea::GetCollectionAreaAbs() const
     123{
     124    return TMath::Pi()*fMcAreaRadius*fMcAreaRadius;
    114125}
    115126

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

    r8673 r8911  
    5858    Double_t GetEntries() const { return fHistAll.Integral(); }
    5959    Double_t GetCollectionAreaEff() const { return fHEnergy.Integral(); }
    60     Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; }
     60    Double_t GetCollectionAreaAbs() const;
    6161
    6262    void SetMcAreaRadius(Float_t x) { fMcAreaRadius = x; }

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

    r8047 r8911  
    1818!   Author(s): Thomas Bretz  12/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2005
     20!   Copyright: MAGIC Software Development, 2000-2008
    2121!
    2222!
    … …  
    3535//////////////////////////////////////////////////////////////////////////////
    3636#include "MHThreshold.h"
     37
     38#include <TMath.h>
    3739
    3840#include <TF1.h>

  • trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.cc

    r5691 r8911  
    5555    fName  = name  ? name  : "MHexagonFreqSpace";
    5656    fTitle = title ? title : "Storage container for fourier space camera";
     57}
     58
     59Double_t MHexagonFreqSpace::GetAbs(UInt_t idx) const
     60{
     61    return TMath::Hypot(fDataRe[idx], fDataIm[idx]);
    5762}
    5863

  • trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.h

    r5691 r8911  
    2626    void Set(const MArrayD &re, const MArrayD &im);
    2727
    28     Double_t GetAbs(UInt_t idx) const { return TMath::Hypot(fDataRe[idx], fDataIm[idx]); }
     28    Double_t GetAbs(UInt_t idx) const;
    2929
    3030    static MGeomCam *NewGeomCam(UShort_t num);

  • trunk/MagicSoft/Mars/mimage/MHillas.cc

    r7438 r8911  
    2121!   Author(s): Wolfgang Wittek  6/2002 <mailto:wittek@mppmu.mpg.de>
    2222!
    23 !   Copyright: MAGIC Software Development, 2000-2003
     23!   Copyright: MAGIC Software Development, 2000-2008
    2424!
    2525!
    … …  
    109109// --------------------------------------------------------------------------
    110110//
     111//  return distance to center (coordinate origin)
     112//
     113Double_t MHillas::GetDist0() const
     114{
     115    return TMath::Hypot(fMeanX, fMeanY);
     116}
     117
     118// --------------------------------------------------------------------------
     119//
    111120//  return the mean position as TVector2(meanx, meany)
    112121//

  • trunk/MagicSoft/Mars/mimage/MHillas.h

    r7224 r8911  
    4343    Float_t  GetMeanX() const  { return fMeanX; }
    4444    Float_t  GetMeanY() const  { return fMeanY; }
    45     Double_t GetDist0() const  { return TMath::Hypot(fMeanX, fMeanY); } // return distance to center
     45    Double_t GetDist0() const;
    4646
    4747    TVector2 GetMean() const;

  • trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc

    r8625 r8911  
    1919!   Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
    2020!
    21 !   Copyright: MAGIC Software Development, 2000-2005
     21!   Copyright: MAGIC Software Development, 2000-2008
    2222!
    2323!
    … …  
    9595    fTime      =  0;
    9696    fTimeRms   = -1;
     97}
     98
     99// --------------------------------------------------------------------------
     100//
     101//  return TMath::Hypot(fCenterX, fCenterY);
     102//
     103Float_t MMuonSearchPar::GetDist() const
     104{
     105    return TMath::Hypot(fCenterX, fCenterY);
    97106}
    98107

  • trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h

    r8625 r8911  
    4242    Float_t GetCenterX()   const { return fCenterX; }
    4343    Float_t GetCenterY()   const { return fCenterY; }
    44     Float_t GetDist()      const { return TMath::Hypot(fCenterX, fCenterY); }
     44    Float_t GetDist() const;
    4545    Float_t GetTime()      const { return fTime; }
    4646    Float_t GetTimeRms()   const { return fTimeRms; }

  • trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc

    r8497 r8911  
    2020!   Author(s): Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
    2121!
    22 !   Copyright: MAGIC Software Development, 2000-2004
     22!   Copyright: MAGIC Software Development, 2000-2008
    2323!
    2424!
    … …  
    7373{
    7474    Clear();
     75}
     76
     77// ------------------------------------------------------------------------
     78//
     79//  return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents)   : 0;
     80//
     81Float_t MPedestalPix::GetPedestalError() const
     82{
     83    return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents)   : 0;
     84}
     85
     86// ------------------------------------------------------------------------
     87//
     88// return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents*2) : 0;
     89//
     90Float_t MPedestalPix::GetPedestalRmsError() const
     91{
     92    return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents*2) : 0;
    7593}
    7694

  • trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h

    r5464 r8911  
    3535    Float_t GetPedestalRms() const { return fPedestalRms; }
    3636    Float_t GetPedestalABoffset() const { return fPedestalABoffset; }
    37     Float_t GetPedestalError()    const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents)   : 0; }
    38     Float_t GetPedestalRmsError() const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents*2) : 0; }
     37
     38    Float_t GetPedestalError() const;
     39    Float_t GetPedestalRmsError() const;
    3940
    4041    UInt_t  GetNumEvents() const { return fNumEvents; }

  • trunk/MagicSoft/Mars/mpointing/MPointingDev.cc

    r7203 r8911  
    1818!   Author(s): Thomas Bretz, 07/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2005
     20!   Copyright: MAGIC Software Development, 2000-2008
    2121!
    2222!
    … …  
    3636using namespace std;
    3737
     38Double_t MPointingDev::GetDevZdRad() const
     39{
     40    return fDevZd*TMath::DegToRad();
     41}
     42
     43Double_t MPointingDev::GetDevAzRad() const
     44{
     45    return fDevAz*TMath::DegToRad();
     46}
Note: See TracChangeset for help on using the changeset viewer.