Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc (revision 8910) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc (revision 8911) @@ -18,5 +18,5 @@ ! Author(s): Markus Gaug 02/2004 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -432,19 +432,19 @@ { case kEPoisson4: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6); + fSinglePheFit = new TF1("SinglePheFit",&PoissonKto4,rmin,rmax,6); rmin += 6.5; break; case kEPoisson5: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6); + fSinglePheFit = new TF1("SinglePheFit",&PoissonKto5,rmin,rmax,6); rmin = 0.; break; case kEPoisson6: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6); + fSinglePheFit = new TF1("SinglePheFit",&PoissonKto6,rmin,rmax,6); break; case kEPolya: - fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8); + fSinglePheFit = new TF1("SinglePheFit",&Polya,rmin,rmax,8); break; case kEMichele: - fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9); + fSinglePheFit = new TF1("SinglePheFit",&FitFuncMichele,rmin,rmax,9); break; default: @@ -964,12 +964,341 @@ } - - - - - - - - - - +Double_t MHCalibrationChargeBlindPix::FitFuncMichele(Double_t *x, Double_t *par) +{ + Double_t lambda1cat = par[0]; + Double_t lambda1dyn = par[1]; + Double_t mu0 = par[2]; + Double_t mu1cat = par[3]; + Double_t mu1dyn = par[4]; + Double_t sigma0 = par[5]; + Double_t sigma1cat = par[6]; + Double_t sigma1dyn = par[7]; + + Double_t sumcat = 0.; + Double_t sumdyn = 0.; + Double_t arg = 0.; + + if (lambda1cat < lambda1dyn) + return FLT_MAX; + + if (mu1cat < mu0) + return FLT_MAX; + + if (mu1dyn < mu0) + return FLT_MAX; + + if (mu1cat < mu1dyn) + return FLT_MAX; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1cat < sigma0) + return FLT_MAX; + + if (sigma1dyn < sigma0) + return FLT_MAX; + + Double_t mu2cat = (2.*mu1cat)-mu0; + Double_t mu2dyn = (2.*mu1dyn)-mu0; + Double_t mu3cat = (3.*mu1cat)-(2.*mu0); + Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0); + + Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); + Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); + Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0)); + Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0)); + + Double_t lambda2cat = lambda1cat*lambda1cat; + Double_t lambda2dyn = lambda1dyn*lambda1dyn; + Double_t lambda3cat = lambda2cat*lambda1cat; + Double_t lambda3dyn = lambda2dyn*lambda1dyn; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sumcat = TMath::Exp(-0.5*arg*arg)/sigma0; + sumdyn = sumcat; + + // k=1cat: + arg = (x[0] - mu1cat)/sigma1cat; + sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat; + // k=1dyn: + arg = (x[0] - mu1dyn)/sigma1dyn; + sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn; + + // k=2cat: + arg = (x[0] - mu2cat)/sigma2cat; + sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat; + // k=2dyn: + arg = (x[0] - mu2dyn)/sigma2dyn; + sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn; + + + // k=3cat: + arg = (x[0] - mu3cat)/sigma3cat; + sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat; + // k=3dyn: + arg = (x[0] - mu3dyn)/sigma3dyn; + sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; + + sumcat = TMath::Exp(-1.*lambda1cat)*sumcat; + sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn; + + return par[8]*(sumcat+sumdyn)/2.; +} + +Double_t MHCalibrationChargeBlindPix::PoissonKto4(Double_t *x, Double_t *par) +{ + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return FLT_MAX; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1 < sigma0) + return FLT_MAX; + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + return TMath::Exp(-1.*lambda)*par[5]*sum; +} + +Double_t MHCalibrationChargeBlindPix::PoissonKto5(Double_t *x, Double_t *par) +{ + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return FLT_MAX; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1 < sigma0) + return FLT_MAX; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + return TMath::Exp(-1.*lambda)*par[5]*sum; +} + +Double_t MHCalibrationChargeBlindPix::PoissonKto6(Double_t *x, Double_t *par) +{ + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return FLT_MAX; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma0 < 0.0001) + return FLT_MAX; + + if (sigma1 < sigma0) + return FLT_MAX; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + Double_t mu6 = (6.*mu1)-(5.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + Double_t lambda6 = lambda5*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + // k=6: + arg = (x[0] - mu6)/sigma6; + sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; + + return TMath::Exp(-1.*lambda)*par[5]*sum; +} + +Double_t MHCalibrationChargeBlindPix::Polya(Double_t *x, Double_t *par) +{ + const Double_t QEcat = 0.247; // mean quantum efficiency + const Double_t sqrt2 = 1.4142135623731; + const Double_t sqrt3 = 1.7320508075689; + const Double_t sqrt4 = 2.; + + const Double_t lambda = par[0]; // mean number of photons + + const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution + const Double_t delta1 = par[2]; // amplification first dynode + const Double_t delta2 = par[3]; // amplification subsequent dynodes + + const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges + + const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain + const Double_t A = 1. + excessPoisson - QEcat + + 1./delta1 + + 1./delta1/delta2 + + 1./delta1/delta2/delta2; // variance contributions from PMT and QE + + const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion + + const Double_t mu0 = par[7]; // pedestal + const Double_t mu1 = totAmpl; // single phe position + const Double_t mu2 = 2*totAmpl; // double phe position + const Double_t mu3 = 3*totAmpl; // triple phe position + const Double_t mu4 = 4*totAmpl; // quadruple phe position + + const Double_t sigma0 = par[5]; + const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A); + const Double_t sigma2 = sqrt2*sigma1; + const Double_t sigma3 = sqrt3*sigma1; + const Double_t sigma4 = sqrt4*sigma1; + + const Double_t lambda2 = lambda*lambda; + const Double_t lambda3 = lambda2*lambda; + const Double_t lambda4 = lambda3*lambda; + + //-- calculate the area---- + Double_t arg = (x[0] - mu0)/sigma0; + Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + return TMath::Exp(-1.*lambda)*par[6]*sum; +} Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h (revision 8910) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h (revision 8911) @@ -5,4 +5,6 @@ #include "MHCalibrationPix.h" #endif + +#include #ifndef ROOT_TF1 @@ -59,4 +61,10 @@ Byte_t fFlags; // Bit-field for the flags enum { kSinglePheFitOK, kPedestalFitOK }; // Possible bits to be set + + static Double_t FitFuncMichele(Double_t *x, Double_t *par); + static Double_t PoissonKto4(Double_t *x, Double_t *par); + static Double_t PoissonKto5(Double_t *x, Double_t *par); + static Double_t PoissonKto6(Double_t *x, Double_t *par); + static Double_t Polya(Double_t *x, Double_t *par); public: @@ -137,357 +145,4 @@ const Double_t sigma0, const Double_t sigma1); -private: - - inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par) - { - - Double_t lambda1cat = par[0]; - Double_t lambda1dyn = par[1]; - Double_t mu0 = par[2]; - Double_t mu1cat = par[3]; - Double_t mu1dyn = par[4]; - Double_t sigma0 = par[5]; - Double_t sigma1cat = par[6]; - Double_t sigma1dyn = par[7]; - - Double_t sumcat = 0.; - Double_t sumdyn = 0.; - Double_t arg = 0.; - - if (lambda1cat < lambda1dyn) - return FLT_MAX; - - if (mu1cat < mu0) - return FLT_MAX; - - if (mu1dyn < mu0) - return FLT_MAX; - - if (mu1cat < mu1dyn) - return FLT_MAX; - - if (sigma0 < 0.0001) - return FLT_MAX; - - if (sigma1cat < sigma0) - return FLT_MAX; - - if (sigma1dyn < sigma0) - return FLT_MAX; - - Double_t mu2cat = (2.*mu1cat)-mu0; - Double_t mu2dyn = (2.*mu1dyn)-mu0; - Double_t mu3cat = (3.*mu1cat)-(2.*mu0); - Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0); - - Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); - Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); - Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0)); - Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0)); - - Double_t lambda2cat = lambda1cat*lambda1cat; - Double_t lambda2dyn = lambda1dyn*lambda1dyn; - Double_t lambda3cat = lambda2cat*lambda1cat; - Double_t lambda3dyn = lambda2dyn*lambda1dyn; - - // k=0: - arg = (x[0] - mu0)/sigma0; - sumcat = TMath::Exp(-0.5*arg*arg)/sigma0; - sumdyn = sumcat; - - // k=1cat: - arg = (x[0] - mu1cat)/sigma1cat; - sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat; - // k=1dyn: - arg = (x[0] - mu1dyn)/sigma1dyn; - sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn; - - // k=2cat: - arg = (x[0] - mu2cat)/sigma2cat; - sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat; - // k=2dyn: - arg = (x[0] - mu2dyn)/sigma2dyn; - sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn; - - - // k=3cat: - arg = (x[0] - mu3cat)/sigma3cat; - sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat; - // k=3dyn: - arg = (x[0] - mu3dyn)/sigma3dyn; - sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; - - sumcat = TMath::Exp(-1.*lambda1cat)*sumcat; - sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn; - - return par[8]*(sumcat+sumdyn)/2.; - - } - - inline static Double_t fPoissonKto4(Double_t *x, Double_t *par) - { - - Double_t lambda = par[0]; - - Double_t sum = 0.; - Double_t arg = 0.; - - Double_t mu0 = par[1]; - Double_t mu1 = par[2]; - - if (mu1 < mu0) - return FLT_MAX; - - Double_t sigma0 = par[3]; - Double_t sigma1 = par[4]; - - if (sigma0 < 0.0001) - return FLT_MAX; - - if (sigma1 < sigma0) - return FLT_MAX; - - Double_t mu2 = (2.*mu1)-mu0; - Double_t mu3 = (3.*mu1)-(2.*mu0); - Double_t mu4 = (4.*mu1)-(3.*mu0); - - Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); - Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); - Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); - - Double_t lambda2 = lambda*lambda; - Double_t lambda3 = lambda2*lambda; - Double_t lambda4 = lambda3*lambda; - - // k=0: - arg = (x[0] - mu0)/sigma0; - sum = TMath::Exp(-0.5*arg*arg)/sigma0; - - // k=1: - arg = (x[0] - mu1)/sigma1; - sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; - - // k=2: - arg = (x[0] - mu2)/sigma2; - sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; - - // k=3: - arg = (x[0] - mu3)/sigma3; - sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; - - // k=4: - arg = (x[0] - mu4)/sigma4; - sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; - - return TMath::Exp(-1.*lambda)*par[5]*sum; - - } - - - inline static Double_t fPoissonKto5(Double_t *x, Double_t *par) - { - - Double_t lambda = par[0]; - - Double_t sum = 0.; - Double_t arg = 0.; - - Double_t mu0 = par[1]; - Double_t mu1 = par[2]; - - if (mu1 < mu0) - return FLT_MAX; - - Double_t sigma0 = par[3]; - Double_t sigma1 = par[4]; - - if (sigma0 < 0.0001) - return FLT_MAX; - - if (sigma1 < sigma0) - return FLT_MAX; - - - Double_t mu2 = (2.*mu1)-mu0; - Double_t mu3 = (3.*mu1)-(2.*mu0); - Double_t mu4 = (4.*mu1)-(3.*mu0); - Double_t mu5 = (5.*mu1)-(4.*mu0); - - Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); - Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); - Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); - Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); - - Double_t lambda2 = lambda*lambda; - Double_t lambda3 = lambda2*lambda; - Double_t lambda4 = lambda3*lambda; - Double_t lambda5 = lambda4*lambda; - - // k=0: - arg = (x[0] - mu0)/sigma0; - sum = TMath::Exp(-0.5*arg*arg)/sigma0; - - // k=1: - arg = (x[0] - mu1)/sigma1; - sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; - - // k=2: - arg = (x[0] - mu2)/sigma2; - sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; - - // k=3: - arg = (x[0] - mu3)/sigma3; - sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; - - // k=4: - arg = (x[0] - mu4)/sigma4; - sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; - - // k=5: - arg = (x[0] - mu5)/sigma5; - sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; - - return TMath::Exp(-1.*lambda)*par[5]*sum; - - } - - - inline static Double_t fPoissonKto6(Double_t *x, Double_t *par) - { - - Double_t lambda = par[0]; - - Double_t sum = 0.; - Double_t arg = 0.; - - Double_t mu0 = par[1]; - Double_t mu1 = par[2]; - - if (mu1 < mu0) - return FLT_MAX; - - Double_t sigma0 = par[3]; - Double_t sigma1 = par[4]; - - if (sigma0 < 0.0001) - return FLT_MAX; - - if (sigma1 < sigma0) - return FLT_MAX; - - - Double_t mu2 = (2.*mu1)-mu0; - Double_t mu3 = (3.*mu1)-(2.*mu0); - Double_t mu4 = (4.*mu1)-(3.*mu0); - Double_t mu5 = (5.*mu1)-(4.*mu0); - Double_t mu6 = (6.*mu1)-(5.*mu0); - - Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); - Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); - Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); - Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); - Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); - - Double_t lambda2 = lambda*lambda; - Double_t lambda3 = lambda2*lambda; - Double_t lambda4 = lambda3*lambda; - Double_t lambda5 = lambda4*lambda; - Double_t lambda6 = lambda5*lambda; - - // k=0: - arg = (x[0] - mu0)/sigma0; - sum = TMath::Exp(-0.5*arg*arg)/sigma0; - - // k=1: - arg = (x[0] - mu1)/sigma1; - sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; - - // k=2: - arg = (x[0] - mu2)/sigma2; - sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; - - // k=3: - arg = (x[0] - mu3)/sigma3; - sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; - - // k=4: - arg = (x[0] - mu4)/sigma4; - sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; - - // k=5: - arg = (x[0] - mu5)/sigma5; - sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; - - // k=6: - arg = (x[0] - mu6)/sigma6; - sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; - - return TMath::Exp(-1.*lambda)*par[5]*sum; - - } - - inline static Double_t fPolya(Double_t *x, Double_t *par) - { - - const Double_t QEcat = 0.247; // mean quantum efficiency - const Double_t sqrt2 = 1.4142135623731; - const Double_t sqrt3 = 1.7320508075689; - const Double_t sqrt4 = 2.; - - const Double_t lambda = par[0]; // mean number of photons - - const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution - const Double_t delta1 = par[2]; // amplification first dynode - const Double_t delta2 = par[3]; // amplification subsequent dynodes - - const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges - - const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain - const Double_t A = 1. + excessPoisson - QEcat - + 1./delta1 - + 1./delta1/delta2 - + 1./delta1/delta2/delta2; // variance contributions from PMT and QE - - const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion - - const Double_t mu0 = par[7]; // pedestal - const Double_t mu1 = totAmpl; // single phe position - const Double_t mu2 = 2*totAmpl; // double phe position - const Double_t mu3 = 3*totAmpl; // triple phe position - const Double_t mu4 = 4*totAmpl; // quadruple phe position - - const Double_t sigma0 = par[5]; - const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A); - const Double_t sigma2 = sqrt2*sigma1; - const Double_t sigma3 = sqrt3*sigma1; - const Double_t sigma4 = sqrt4*sigma1; - - const Double_t lambda2 = lambda*lambda; - const Double_t lambda3 = lambda2*lambda; - const Double_t lambda4 = lambda3*lambda; - - //-- calculate the area---- - Double_t arg = (x[0] - mu0)/sigma0; - Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0; - - // k=1: - arg = (x[0] - mu1)/sigma1; - sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; - - // k=2: - arg = (x[0] - mu2)/sigma2; - sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; - - // k=3: - arg = (x[0] - mu3)/sigma3; - sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; - - // k=4: - arg = (x[0] - mu4)/sigma4; - sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; - - return TMath::Exp(-1.*lambda)*par[6]*sum; - } - ClassDef(MHCalibrationChargeBlindPix, 1) // Histogram class for Charge Blind Pixel Calibration }; Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 8910) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 8911) @@ -19,5 +19,5 @@ ! Author(s): Harald Kornmayer 1/2001 ! -! Copyright: MAGIC Software Development, 2000-2006 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -35,4 +35,6 @@ ////////////////////////////////////////////////////////////////////////////// #include "MHCollectionArea.h" + +#include #include @@ -112,4 +114,13 @@ fHistAll.Sumw2(); fHEnergy.Sumw2(); +} + +// -------------------------------------------------------------------------- +// +// Return the Area defined by fMcAreaRadius +// +Double_t MHCollectionArea::GetCollectionAreaAbs() const +{ + return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; } Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 8910) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 8911) @@ -58,5 +58,5 @@ Double_t GetEntries() const { return fHistAll.Integral(); } Double_t GetCollectionAreaEff() const { return fHEnergy.Integral(); } - Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; } + Double_t GetCollectionAreaAbs() const; void SetMcAreaRadius(Float_t x) { fMcAreaRadius = x; } Index: trunk/MagicSoft/Mars/mhflux/MHThreshold.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHThreshold.cc (revision 8910) +++ trunk/MagicSoft/Mars/mhflux/MHThreshold.cc (revision 8911) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz 12/2005 ! -! Copyright: MAGIC Software Development, 2000-2005 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -35,4 +35,6 @@ ////////////////////////////////////////////////////////////////////////////// #include "MHThreshold.h" + +#include #include Index: trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.cc =================================================================== --- trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.cc (revision 8910) +++ trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.cc (revision 8911) @@ -55,4 +55,9 @@ fName = name ? name : "MHexagonFreqSpace"; fTitle = title ? title : "Storage container for fourier space camera"; +} + +Double_t MHexagonFreqSpace::GetAbs(UInt_t idx) const +{ + return TMath::Hypot(fDataRe[idx], fDataIm[idx]); } Index: trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.h =================================================================== --- trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.h (revision 8910) +++ trunk/MagicSoft/Mars/mhft/MHexagonFreqSpace.h (revision 8911) @@ -26,5 +26,5 @@ void Set(const MArrayD &re, const MArrayD &im); - Double_t GetAbs(UInt_t idx) const { return TMath::Hypot(fDataRe[idx], fDataIm[idx]); } + Double_t GetAbs(UInt_t idx) const; static MGeomCam *NewGeomCam(UShort_t num); Index: trunk/MagicSoft/Mars/mimage/MHillas.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MHillas.cc (revision 8910) +++ trunk/MagicSoft/Mars/mimage/MHillas.cc (revision 8911) @@ -21,5 +21,5 @@ ! Author(s): Wolfgang Wittek 6/2002 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -109,4 +109,13 @@ // -------------------------------------------------------------------------- // +// return distance to center (coordinate origin) +// +Double_t MHillas::GetDist0() const +{ + return TMath::Hypot(fMeanX, fMeanY); +} + +// -------------------------------------------------------------------------- +// // return the mean position as TVector2(meanx, meany) // Index: trunk/MagicSoft/Mars/mimage/MHillas.h =================================================================== --- trunk/MagicSoft/Mars/mimage/MHillas.h (revision 8910) +++ trunk/MagicSoft/Mars/mimage/MHillas.h (revision 8911) @@ -43,5 +43,5 @@ Float_t GetMeanX() const { return fMeanX; } Float_t GetMeanY() const { return fMeanY; } - Double_t GetDist0() const { return TMath::Hypot(fMeanX, fMeanY); } // return distance to center + Double_t GetDist0() const; TVector2 GetMean() const; Index: trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc (revision 8910) +++ trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc (revision 8911) @@ -19,5 +19,5 @@ ! Author(s): Markus Meyer 10/2004 ! -! Copyright: MAGIC Software Development, 2000-2005 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -95,4 +95,13 @@ fTime = 0; fTimeRms = -1; +} + +// -------------------------------------------------------------------------- +// +// return TMath::Hypot(fCenterX, fCenterY); +// +Float_t MMuonSearchPar::GetDist() const +{ + return TMath::Hypot(fCenterX, fCenterY); } Index: trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h (revision 8910) +++ trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h (revision 8911) @@ -42,5 +42,5 @@ Float_t GetCenterX() const { return fCenterX; } Float_t GetCenterY() const { return fCenterY; } - Float_t GetDist() const { return TMath::Hypot(fCenterX, fCenterY); } + Float_t GetDist() const; Float_t GetTime() const { return fTime; } Float_t GetTimeRms() const { return fTimeRms; } Index: trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc (revision 8910) +++ trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc (revision 8911) @@ -20,5 +20,5 @@ ! Author(s): Florian Goebel 06/2004 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -73,4 +73,22 @@ { Clear(); +} + +// ------------------------------------------------------------------------ +// +// return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents) : 0; +// +Float_t MPedestalPix::GetPedestalError() const +{ + return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents) : 0; +} + +// ------------------------------------------------------------------------ +// +// return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents*2) : 0; +// +Float_t MPedestalPix::GetPedestalRmsError() const +{ + return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents*2) : 0; } Index: trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h (revision 8910) +++ trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h (revision 8911) @@ -35,6 +35,7 @@ Float_t GetPedestalRms() const { return fPedestalRms; } Float_t GetPedestalABoffset() const { return fPedestalABoffset; } - Float_t GetPedestalError() const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents) : 0; } - Float_t GetPedestalRmsError() const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents*2) : 0; } + + Float_t GetPedestalError() const; + Float_t GetPedestalRmsError() const; UInt_t GetNumEvents() const { return fNumEvents; } Index: trunk/MagicSoft/Mars/mpointing/MPointingDev.cc =================================================================== --- trunk/MagicSoft/Mars/mpointing/MPointingDev.cc (revision 8910) +++ trunk/MagicSoft/Mars/mpointing/MPointingDev.cc (revision 8911) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 07/2005 ! -! Copyright: MAGIC Software Development, 2000-2005 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -36,2 +36,11 @@ using namespace std; +Double_t MPointingDev::GetDevZdRad() const +{ + return fDevZd*TMath::DegToRad(); +} + +Double_t MPointingDev::GetDevAzRad() const +{ + return fDevAz*TMath::DegToRad(); +}