Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2912) +++ trunk/MagicSoft/Mars/Changelog (revision 2913) @@ -14,4 +14,6 @@ - TSpline5 now on stack. + * mcalib/MHCalibrationBlindPixel.[h,cc]: + - force mu_0 in Blind Pixel Fit to be around 0 in fKPoisson4 2004/01/26: Thomas Bretz Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 2912) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 2913) @@ -347,12 +347,40 @@ const Stat_t entries = fHBlindPixelCharge->Integral("width"); const Double_t lambda_guess = 0.5; - const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin()); + const Double_t maximum_bin = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin()); + const Double_t norm = entries/gkSq2Pi; + + // + // Initialize the fit function + // + switch (fFitFunc) + { + case kEPoisson4: + // fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,5); + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6); + break; + case kEPoisson5: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6); + break; + case kEPoisson6: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6); + break; + case kEPolya: + break; + + case kEMichele: + break; + + default: + *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; + return kFALSE; + break; + } + + const Double_t mu_0_guess = maximum_bin; const Double_t si_0_guess = 20.; const Double_t mu_1_guess = mu_0_guess + 50.; const Double_t si_1_guess = si_0_guess + si_0_guess; - const Double_t norm = entries/gkSq2Pi; - - // - // Initialize + // + // Initialize boundaries and start parameters // switch (fFitFunc) @@ -360,5 +388,15 @@ case kEPoisson4: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6); + fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); + fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area"); + fSinglePheFit->SetParLimits(0,0.,1.); + fSinglePheFit->SetParLimits(1,-2.,2.); + fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min))); + fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); + fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); + fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5); + break; + case kEPoisson5: + case kEPoisson6: fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area"); @@ -369,31 +407,4 @@ fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); - break; - - case kEPoisson5: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6); - fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); - fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area"); - fSinglePheFit->SetParLimits(0,0.,1.); - fSinglePheFit->SetParLimits(1,min,(max-min)/1.5); - fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min))); - fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); - fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); - fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); - break; - - case kEPoisson6: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6); - fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); - fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area"); - fSinglePheFit->SetParLimits(0,0.,1.); - fSinglePheFit->SetParLimits(1,min,(max-min)/1.5); - fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min))); - fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); - fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); - fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); - break; - - case kEPoisson7: break; Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 2912) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 2913) @@ -132,52 +132,13 @@ { + const Double_t qe = 0.247; + Double_t lambda = par[0]; + + Double_t b1 = par[1]; + Double_t b2 = par[2]; Double_t sum = 0.; Double_t arg = 0.; - - Double_t mu0 = par[1]; - Double_t mu1 = par[2]; - - if (mu1 < mu0) - return fNoWay; - - Double_t sigma0 = par[3]; - Double_t sigma1 = par[4]; - - if (sigma1 < sigma0) - return fNoWay; - - 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; @@ -193,4 +154,5 @@ Double_t arg = 0.; + // Double_t mu0 = 0.; Double_t mu0 = par[1]; Double_t mu1 = par[2];