Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2919) +++ trunk/MagicSoft/Mars/Changelog (revision 2920) @@ -11,18 +11,25 @@ in arrays: gsNCells, gsNTrigPixels, gsNPixInCell, gsNLutInCell, gsNPixInLut, fNumPixCell. - - Added method MMcTriggerLvl2::GetCellCompactPixel(int cell, MGeomCam *fCam) + - Added method + MMcTriggerLvl2::GetCellCompactPixel(int cell, MGeomCam *fCam) which computes compact pixels into a given L2T macrocell. - - Added method MMcTriggerLvl2::CalcBiggerCellPseudoSize() which - computes fCellPseudoSize, the maximum Pseudo Size into L2T macrocells - - Added method MMcTriggerLvl2::GetCellPseudoSize() const which - returns fCellPseudoSize - - Added method MMcTriggerLvl2::IsPixelInCell(Int_t pixel, Int_t cell), + - Added method + MMcTriggerLvl2::CalcBiggerCellPseudoSize() + which computes fCellPseudoSize, the maximum Pseudo Size into L2T + macrocells + - Added method + MMcTriggerLvl2::GetCellPseudoSize() const + which returns fCellPseudoSize + - Added method + MMcTriggerLvl2::IsPixelInCell(Int_t pixel, Int_t cell), which controls whether a pixel belongs to a given L2T cell. - - Added method MMcTriggerLvl2::GetMaxCell() const which returns fMaxCell, the cell with - the maximum fCellPseudoSize. - - - - 2004/01/26: Markus Gaug + - Added method + MMcTriggerLvl2::GetMaxCell() const + which returns fMaxCell, the cell with the maximum + fCellPseudoSize. + + + + 2004/01/26: Markus Gaug, Michele Doro * manalysis/MArrivalTime.[h,cc], manalysis/MArrivalTimeCalc.[h,cc]: @@ -33,5 +40,7 @@ * mcalib/MHCalibrationBlindPixel.[h,cc]: - - force mu_0 in Blind Pixel Fit to be around 0 in fKPoisson4 + - force mu_{0} in Blind Pixel Fit to be around 0 in fKPoisson4 + - implement combined Polya fit and Michele's back-scattered electron + fit Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 2919) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 2920) @@ -260,5 +260,5 @@ { - gStyle->SetOptFit(0); + gStyle->SetOptFit(1); gStyle->SetOptStat(111111); @@ -366,6 +366,6 @@ break; case kEPolya: - break; - + fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8); + break; case kEMichele: break; @@ -381,4 +381,17 @@ const Double_t mu_1_guess = mu_0_guess + 50.; const Double_t si_1_guess = si_0_guess + si_0_guess; + // Michele + const Double_t lambda_1cat_guess = 0.5; + const Double_t lambda_1dyn_guess = 0.5; + const Double_t mu_1cat_guess = mu_0_guess + 50.; + const Double_t mu_1dyn_guess = mu_0_guess + 20.; + const Double_t si_1cat_guess = si_0_guess + si_0_guess; + const Double_t si_1dyn_guess = si_0_guess; + // Polya + const Double_t excessPoisson_guess = 0.5; + const Double_t delta1_guess = 8.; + const Double_t delta2_guess = 5.; + const Double_t electronicAmp_guess = 0.0025; + // // Initialize boundaries and start parameters @@ -389,5 +402,5 @@ case kEPoisson4: 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->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area"); fSinglePheFit->SetParLimits(0,0.,1.); fSinglePheFit->SetParLimits(1,-2.,2.); @@ -400,5 +413,5 @@ 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"); + 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); @@ -410,7 +423,26 @@ case kEPolya: + fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess, + delta1_guess,delta2_guess, + electronicAmp_guess, + si_0_guess, + norm, mu_0_guess); + fSinglePheFit->SetParNames("#lambda","b_{tot}", + "#delta_{1}","#delta_{2}", + "amp_{e}","#sigma_{0}", + "Area", "#mu_{0}"); + fSinglePheFit->SetParLimits(0,0.,1.); + fSinglePheFit->SetParLimits(1,0.,1.); + fSinglePheFit->SetParLimits(2,6.,12.); + fSinglePheFit->SetParLimits(3,3.,8.); + fSinglePheFit->SetParLimits(4,0.,0.005); + fSinglePheFit->SetParLimits(5,min,(max-min)/1.5); + fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1); + fSinglePheFit->SetParLimits(7,-35.,15.); break; case kEMichele: + + break; @@ -450,4 +482,16 @@ fSigma1Err = fSinglePheFit->GetParError(4); break; + case kEPolya: + fLambda = fSinglePheFit->GetParameter(0); + fMu0 = fSinglePheFit->GetParameter(7); + fMu1 = 0.; + fSigma0 = fSinglePheFit->GetParameter(5); + fSigma1 = 0.; + + fLambdaErr = fSinglePheFit->GetParError(0); + fMu0Err = fSinglePheFit->GetParError(7); + fMu1Err = 0.; + fSigma0Err = fSinglePheFit->GetParError(5); + fSigma1Err = 0.; default: break; @@ -479,5 +523,5 @@ // Perform the cross-check fitting only the pedestal: - fSinglePhePedFit = new TF1("GausPed","gaus",rmin,fMu0); + fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.); fHBlindPixelCharge->Fit(fSinglePhePedFit,opt); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 2919) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 2920) @@ -71,31 +71,31 @@ //Getters - const Double_t GetLambda() const { return fLambda; } - const Double_t GetLambdaCheck() const { return fLambdaCheck; } - const Double_t GetMu0() const { return fMu0; } - const Double_t GetMu1() const { return fMu1; } - const Double_t GetSigma0() const { return fSigma0; } - const Double_t GetSigma1() const { return fSigma1; } - - const Double_t GetLambdaErr() const { return fLambdaErr; } - const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; } - const Double_t GetMu0Err() const { return fMu0Err; } - const Double_t GetMu1Err() const { return fMu1Err; } - const Double_t GetSigma0Err() const { return fSigma0Err; } - const Double_t GetSigma1Err() const { return fSigma1Err; } - - const Double_t GetChiSquare() const { return fChisquare; } - const Double_t GetProb() const { return fProb; } - const Int_t GetNdf() const { return fNdf; } - - const Double_t GetMeanTime() const { return fMeanTime; } - const Double_t GetMeanTimeErr() const { return fMeanTimeErr; } - const Double_t GetSigmaTime() const { return fSigmaTime; } - const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; } - - const Bool_t IsFitOK() { return fFitOK; } - - const TH1I *GetHBlindPixelChargevsN() const { return fHBlindPixelChargevsN; } - const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD; } + const Double_t GetLambda() const { return fLambda; } + const Double_t GetLambdaCheck() const { return fLambdaCheck; } + const Double_t GetMu0() const { return fMu0; } + const Double_t GetMu1() const { return fMu1; } + const Double_t GetSigma0() const { return fSigma0; } + const Double_t GetSigma1() const { return fSigma1; } + + const Double_t GetLambdaErr() const { return fLambdaErr; } + const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; } + const Double_t GetMu0Err() const { return fMu0Err; } + const Double_t GetMu1Err() const { return fMu1Err; } + const Double_t GetSigma0Err() const { return fSigma0Err; } + const Double_t GetSigma1Err() const { return fSigma1Err; } + + const Double_t GetChiSquare() const { return fChisquare; } + const Double_t GetProb() const { return fProb; } + const Int_t GetNdf() const { return fNdf; } + + const Double_t GetMeanTime() const { return fMeanTime; } + const Double_t GetMeanTimeErr() const { return fMeanTimeErr; } + const Double_t GetSigmaTime() const { return fSigmaTime; } + const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; } + + const Bool_t IsFitOK() const { return fFitOK; } + + const TH1I *GetHBlindPixelChargevsN() const { return fHBlindPixelChargevsN; } + const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD; } // Draws @@ -131,17 +131,80 @@ inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par) { - - 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 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.; - return TMath::Exp(-1.*lambda)*par[5]*sum; - + if (mu1cat < mu0) + return fNoWay; + + if (sigma1cat < sigma0) + return fNoWay; + + // if (sigma1cat < sigma1dyn) + // return NoWay; + + //if (mu1cat < mu1dyn) + // return NoWay; + + // if (lambda1cat < lambda1dyn) + // return NoWay; + + 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.; + } @@ -339,4 +402,69 @@ } + + 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(MHCalibrationBlindPixel, 1) // Histograms from the Calibration Blind Pixel