Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 4197) +++ trunk/MagicSoft/Mars/Changelog (revision 4198) @@ -29,4 +29,8 @@ - set MCalibrationChargeBlindPix and MCalibrationChargePINDiode as data members and store them in the output of WriteResult() + + * mcalib/MCalibrationQEPix.cc + - several small bugs fixed concerning calibration with the blind pixel + or the PIN Diode * mjobs/MGCamDisplays.[h,cc] Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc (revision 4197) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc (revision 4198) @@ -227,5 +227,5 @@ const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit = 1.; const Float_t MCalibrationChargeCalc::fgLambdaErrLimit = 0.2; -const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit = 0.2; +const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit = 0.5; const Float_t MCalibrationChargeCalc::fgPheErrLimit = 3.5; const Float_t MCalibrationChargeCalc::fgFFactorErrLimit = 3.; @@ -600,5 +600,5 @@ *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run " << "instead of a calibration run " << endl; - return kFALSE; + return kFALSE; } @@ -1117,15 +1117,18 @@ const Float_t lambdacheck = fBlindPixel->GetLambdaCheck(); - if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) < fLambdaCheckLimit) - { - *fLog << warn << GetDescriptor() << ": Lambda and Lambda-Check differ by more than " - << fLambdaCheckLimit << " in the Blind Pixel " << endl; - return kFALSE; - } - - if (lambdaerr < fLambdaErrLimit) - { - *fLog << warn << GetDescriptor() << ": Error of Fitted Lambda is greater than " - << fLambdaErrLimit << " in Blind Pixel " << endl; + if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit) + { + *fLog << warn << GetDescriptor() + << Form("%s%4.2f%s%4.2f%s%4.2f%s",": Lambda: ",lambda," and Lambda-Check: ", + lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel ") + << endl; + return kFALSE; + } + + if (lambdaerr > fLambdaErrLimit) + { + *fLog << warn << GetDescriptor() + << Form("%s%4.2f%s%4.2f%s",": Error of Fitted Lambda: ",lambdaerr," is greater than ", + fLambdaErrLimit," in Blind Pixel ") << endl; return kFALSE; } @@ -1376,5 +1379,5 @@ continue; } - + if (!fBlindPixel->IsFluxInsidePlexiglassAvailable()) { @@ -1385,5 +1388,5 @@ MBadPixelsPix &bad = (*fBadPixels)[i]; - if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun)) + if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun)) { qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor); @@ -1397,5 +1400,5 @@ / fBlindPixel->GetFluxInsidePlexiglass() / geo.GetA() - * fQECam->GetPlexiglassQE(); + * fQECam->GetPlexiglassQE(); const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar() @@ -1405,5 +1408,9 @@ qepix.SetQEBlindPixel ( qe , fPulserColor ); qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor ); - qepix.UpdateBlindPixelMethod(); + qepix.SetBlindPixelMethodValid( kTRUE , fPulserColor ); + + if (!qepix.UpdateBlindPixelMethod()) + *fLog << warn << GetDescriptor() + << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl; } } @@ -1471,5 +1478,9 @@ qepix.SetQEPINDiode ( qe , fPulserColor ); qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor ); - qepix.UpdateBlindPixelMethod(); + qepix.SetPINDiodeMethodValid( kTRUE , fPulserColor ); + + if (!qepix.UpdatePINDiodeMethod()) + *fLog << warn << GetDescriptor() + << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl; } } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 4197) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 4198) @@ -96,8 +96,8 @@ const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002; -const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 300; -const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -100.5; -const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 5199.5; -const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 200.; +const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 1000; +const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -0.5; +const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 5999.5; +const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 600.; const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.; // -------------------------------------------------------------------------- @@ -146,5 +146,5 @@ SetNumSinglePheLimit(); - SetBinsAfterStripping(30); + SetBinsAfterStripping(50); fHGausHist.SetName("HCalibrationChargeBlindPix"); @@ -216,5 +216,6 @@ fLambdaCheckErr = -999.; - fFitFunc = kEPoisson5; + fFitFunc = kEMichele; + // fFitFunc = kEPoisson4; fNumSinglePhes = 0; @@ -387,6 +388,8 @@ // const Float_t signal = (Float_t)fSignal->GetExtractedSignal(); - if (signal > -1) + if (signal > -0.5) FillHistAndArray(signal); + else + return kTRUE; // @@ -448,8 +451,8 @@ fBlindPix->SetSinglePheFitOK(); else - fBlindPix->SetValid(kFALSE); - - fBlindPix->SetLambda ( fLambda ); - fBlindPix->SetLambdaVar ( fLambdaErr*fLambdaErr ); + fBlindPix->SetValid(IsPedestalFitOK()); + + fBlindPix->SetLambda ( fLambdaCheck ); + fBlindPix->SetLambdaVar ( fLambdaCheckErr*fLambdaCheckErr ); fBlindPix->SetMu0 ( fMu0 ); fBlindPix->SetMu0Err ( fMu0Err ); @@ -462,6 +465,6 @@ fBlindPix->SetProb ( fProb ); - fBlindPix->SetLambdaCheck ( fLambdaCheck ); - fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr ); + fBlindPix->SetLambdaCheck ( fLambda ); + fBlindPix->SetLambdaCheckErr ( fLambdaErr ); return kTRUE; @@ -598,6 +601,6 @@ break; case kEMichele: - break; - + fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,10); + break; default: *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; @@ -614,13 +617,16 @@ const Double_t mu_0_guess = maximum_bin; const Double_t si_0_guess = 40.; - const Double_t mu_1_guess = mu_0_guess + 100.; + const Double_t mu_1_guess = mu_0_guess + 4000.; 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; + const Double_t lambda_1cat_guess = 0.005; + const Double_t lambda_1dyn_guess = 0.002; + // const Double_t mu_1cat_guess = mu_0_guess + 4500.; + // const Double_t mu_1dyn_guess = mu_0_guess + 2500.; + const Double_t mu_1cat_guess = mu_0_guess + 1000.; + const Double_t mu_1dyn_guess = mu_0_guess + 2500.; + const Double_t si_1cat_guess = 500.; + const Double_t si_1dyn_guess = 1000.; + const Double_t offset_guess = 0.5; // Polya const Double_t excessPoisson_guess = 0.5; @@ -689,4 +695,26 @@ break; case kEMichele: + fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess, + mu_0_guess, mu_1cat_guess,mu_1dyn_guess, + si_0_guess, si_1cat_guess,si_1dyn_guess, + norm, offset_guess); + fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}", + "#mu_{0}","#mu_{1cat}","#mu_{1dyn}", + "#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}", + "Area","offset"); + fSinglePheFit->SetParLimits(0,0.,0.2); + fSinglePheFit->SetParLimits(1,0.,0.05); + fSinglePheFit->SetParLimits(2,0.,fSinglePheCut); + // fSinglePheFit->SetParLimits(3,4000.,5500.); + // fSinglePheFit->SetParLimits(4,2500.,5000.); + fSinglePheFit->SetParLimits(3,fSinglePheCut,3000.); + fSinglePheFit->SetParLimits(4,fSinglePheCut,3500.); + fSinglePheFit->SetParLimits(5,0.,fSinglePheCut); + fSinglePheFit->SetParLimits(6,500.,1000.); + fSinglePheFit->SetParLimits(7,300.,1500.); + // fSinglePheFit->SetParLimits(6,300.,2000.); + // fSinglePheFit->SetParLimits(7,100.,2500.); + fSinglePheFit->SetParLimits(8,norm/1.1,norm*1.1); + fSinglePheFit->SetParLimits(9,0.,1.); break; @@ -745,4 +773,17 @@ fSigma0Err = fSinglePheFit->GetParError(5); fSigma1Err = 0.; + case kEMichele: + fLambda = fSinglePheFit->GetParameter(0); + fMu0 = fSinglePheFit->GetParameter(2); + fMu1 = fSinglePheFit->GetParameter(3); + fSigma0 = fSinglePheFit->GetParameter(5); + fSigma1 = fSinglePheFit->GetParameter(6); + + fLambdaErr = fSinglePheFit->GetParError(0); + fMu0Err = fSinglePheFit->GetParError(2); + fMu1Err = fSinglePheFit->GetParError(3); + fSigma0Err = fSinglePheFit->GetParError(5); + fSigma1Err = fSinglePheFit->GetParError(6); + break; default: break; @@ -832,18 +873,17 @@ // Perform the cross-check fitting only the pedestal: const Axis_t rmin = 0.; - const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + // const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Axis_t rmax = fSinglePheCut; FitGaus(opt, rmin, rmax); const Stat_t entries = fHGausHist.Integral("width"); - const Double_t fitarea = fFGausFit->GetParameter(0); - const Double_t pedarea = fitarea * TMath::Sqrt(TMath::TwoPi()) * fFGausFit->GetParameter(2); - + const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut); + fLambdaCheck = TMath::Log(entries/pedarea); - fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0) - + fFGausFit->GetParError(2)/fFGausFit->GetParameter(2); - - - SetPedestalFitOK(); + // estimate the error by the error of the obtained area from the Gauss-function: + fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0); + + SetPedestalFitOK(IsGausFitOK()); return; } @@ -861,9 +901,9 @@ fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (", - (fFitFunc = kEPoisson4) ? "Poisson(k=4))" : - (fFitFunc = kEPoisson5) ? "Poisson(k=5))" : - (fFitFunc = kEPoisson6) ? "Poisson(k=4))" : - (fFitFunc = kEPolya ) ? "Polya(k=4))" : - (fFitFunc = kEMichele ) ? "Michele)" : " none )" )); + (fFitFunc == kEPoisson4) ? "Poisson(k=4))" : + (fFitFunc == kEPoisson5) ? "Poisson(k=5))" : + (fFitFunc == kEPoisson6) ? "Poisson(k=6))" : + (fFitFunc == kEPolya ) ? "Polya(k=4))" : + (fFitFunc == kEMichele ) ? "Michele)" : " none )" )); fFitLegend->SetTextSize(0.05); } @@ -979,4 +1019,9 @@ fFGausFit->SetLineColor(kBlue); fFGausFit->Draw("same"); + TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.); + line->SetBit(kCanDelete); + line->SetLineColor(kBlue); + line->SetLineWidth(3); + line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.); } if (fSinglePheFit) @@ -985,4 +1030,5 @@ fSinglePheFit->Draw("same"); } + pad->cd(2); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h (revision 4197) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h (revision 4198) @@ -146,5 +146,5 @@ // Fits public: - enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; // The possible fit functions + enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; // The possible fit functions private: @@ -175,5 +175,6 @@ Double_t sigma1cat = par[6]; Double_t sigma1dyn = par[7]; - + Double_t offset = par[9]; + Double_t sumcat = 0.; Double_t sumdyn = 0.; @@ -240,5 +241,5 @@ sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn; - return par[8]*(sumcat+sumdyn)/2.; + return par[8]*(sumcat+sumdyn)/2. + offset; } Index: trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 4197) +++ trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 4198) @@ -662,14 +662,15 @@ StripZeros(&fHGausHist,fBinsAfterStripping); + TAxis *axe = fHGausHist.GetXaxis(); // // Get the fitting ranges // - Axis_t rmin = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()) : xmin; - Axis_t rmax = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) : xmax; + Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin; + Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax; // // First guesses for the fit (should be as close to reality as possible, // - const Stat_t entries = fHGausHist.Integral("width"); + const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width"); const Double_t mu_guess = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); const Double_t sigma_guess = fHGausHist.GetRMS(); @@ -687,5 +688,5 @@ fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess); fFGausFit->SetParNames("Area","#mu","#sigma"); - fFGausFit->SetParLimits(0,0.,area_guess*1.5); + fFGausFit->SetParLimits(0,0.,area_guess*25.); fFGausFit->SetParLimits(1,rmin,rmax); fFGausFit->SetParLimits(2,0.,rmax-rmin);