Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 5593) +++ trunk/MagicSoft/Mars/Changelog (revision 5594) @@ -28,4 +28,7 @@ * msignal/MExtractTimeAndChargeSpline.cc - fixed some smaller bugs for special cases + + * mcalib/MCalibrationIntensityChargeCam.[h,cc] + - implemented automatic plotting of Razmik plot and others Index: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc (revision 5593) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc (revision 5594) @@ -41,6 +41,13 @@ #include "MCalibrationIntensityChargeCam.h" #include "MCalibrationChargeCam.h" +#include "MCalibrationChargePix.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" #include +#include +#include +#include ClassImp(MCalibrationIntensityChargeCam); @@ -71,2 +78,246 @@ } +// ------------------------------------------------------------------- +// +// Returns a TGraphErrors with the number of photo-electrons vs. +// the extracted signal of pixel "pixid". +// +TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col) +{ + + const Int_t size = GetSize(); + + TArrayF phe(size); + TArrayF pheerr(size); + TArrayF sig(size); + TArrayF sigerr(size); + + for (Int_t i=0;iGetPulserColor() != col) + continue; + // + // Get the calibration pix from the calibration cam + // + MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid]; + // + // Don't use bad pixels + // + if (!pix.IsFFactorMethodValid()) + continue; + // + phe[i] = pix.GetPheFFactorMethod(); + pheerr[i] = pix.GetPheFFactorMethodErr(); + // + // For the calculation of Q, we have to use the + // converted value! + // + sig [i] = pix.GetConvertedMean(); + sigerr[i] = pix.GetConvertedMeanErr(); + } + + TGraphErrors *gr = new TGraphErrors(size, + sig.GetArray(),phe.GetArray(), + sigerr.GetArray(),pheerr.GetArray()); + gr->SetTitle(Form("%s%3i","Pixel ",pixid)); + gr->GetXaxis()->SetTitle("Q [FADC counts]"); + gr->GetYaxis()->SetTitle("photo-electrons [1]"); + return gr; +} + +// ------------------------------------------------------------------- +// +// Returns a TGraphErrors with the number of photo-electrons vs. +// the extracted signal over all pixels with area index "aidx". +// +// The points represent the means of the pixels values, while the error bars +// the sigma of the pixels values. +// +TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col) +{ + + const Int_t size = GetSize(); + + TArrayF phe(size); + TArrayF pheerr(size); + TArrayF sig(size); + TArrayF sigerr(size); + + for (Int_t i=0;iGetPulserColor() != col) + continue; + + // + // Get the area calibration pix from the calibration cam + // + MCalibrationChargePix &pix = (MCalibrationChargePix&)(cam->GetAverageArea(aidx)); + + phe[i] = pix.GetPheFFactorMethod(); + pheerr[i] = pix.GetPheFFactorMethodErr(); + // + // For the calculation of Q, we have to use the + // converted value! + // + sig [i] = pix.GetConvertedMean(); + sigerr[i] = pix.GetConvertedMeanErr(); + } + + TGraphErrors *gr = new TGraphErrors(size, + sig.GetArray(),phe.GetArray(), + sigerr.GetArray(),pheerr.GetArray()); + gr->SetTitle(Form("%s%3i","Area Index ",aidx)); + gr->GetXaxis()->SetTitle("Q [FADC counts]"); + gr->GetYaxis()->SetTitle("photo-electrons [1]"); + return gr; +} + +// ------------------------------------------------------------------- +// +// Returns a TGraphErrors with the 'Razmik plot' of pixel "pixid". +// The Razmik plot shows the value of 'R' vs. 1/Q where: +// +// sigma^2 F^2 +// R = ------- = ------ +// ^2 +// +// and 1/Q is the inverse (mean) extracted signal +// +TGraphErrors *MCalibrationIntensityChargeCam::GetRazmikPlot( const UInt_t pixid ) +{ + + const Int_t size = GetSize(); + + TArrayF r(size); + TArrayF rerr(size); + TArrayF oneoverq(size); + TArrayF oneoverqerr(size); + + for (Int_t i=0;i^2 ^2 + // + const Float_t mean = pix.GetMean(); + const Float_t meanerr = pix.GetMeanErr(); + const Float_t rsigma = pix.GetRSigma(); + const Float_t rsigmaerr = pix.GetRSigmaErr(); + r[i] = rsigma*rsigma/mean/mean; + const Float_t rrelvar = 4.*rsigmaerr*rsigmaerr/rsigma/rsigma + 4.*meanerr*meanerr/mean/mean; + rerr[i] = rrelvar * r[i] * r[i]; + rerr[i] = rerr[i] <= 0 ? 0. : TMath::Sqrt(rerr[i]); + // + // For the calculation of 1/Q, we have to use the + // converted value! + // + const Float_t q = pix.GetConvertedMean(); + const Float_t qe = pix.GetConvertedMeanErr(); + oneoverq [i] = 1./q; + oneoverqerr[i] = qe / q / q; + } + + TGraphErrors *gr = new TGraphErrors(size, + oneoverq.GetArray(),r.GetArray(), + oneoverqerr.GetArray(),rerr.GetArray()); + gr->SetTitle(Form("%s%3i","Pixel ",pixid)); + gr->GetXaxis()->SetTitle("1/Q [FADC counts^{-1}]"); + gr->GetYaxis()->SetTitle("\sigma_{red}^{2}/Q^{2}"); + return gr; +} + +// ------------------------------------------------------------------- +// +// Returns a 2-dimensional histogram with the fit results of the +// 'Razmik plot' for each pixel of area index "aidx" (see GetRazmikPlot()) +// +// The results of the polynomial fit of grade 1 are: +// +// x-axis: Offset (Parameter 0 of the polynomial) +// y-axis: Slope (Parameter 1 of the polynomial) +// +// The offset is a measure of how well-known the supposed additional contributions +// to the value "reduced sigma" are. Because a photo-multiplier is a linear instrument, +// the excess fluctuations are linear w.r.t. the signal amplitude and can be expressed by +// the proportionality constant F (the "F-Factor"). +// Adding noise from outside (e.g. night sky background) modifies the recorded noise, but +// not the mean extracted signal, due to the AC-coupling. Thus, noise contributions from outside +// (e.g. calculating the pedestal RMS)have to be subtracted from the recorded signal fluctuations +// in order to retrieve the linearity relation: +// +// sigma(signal)^2 / mean(signal)^2 = sigma^2 / ^2 = F^2 / (1) +// +// Any systematic offset in the sigma(signal) will produce an offset in the "Razmik plot"), +// characterized by the Offset of the polynomial fit. Thus, in an ideal case, all pixels have their +// "offset" centered very closely around zero. +// +// The "slope" is the proportionality constant F^2, multiplied with the conversion factor +// phe's to mean signal (because the "Razmik plot" plots the left side of eq. (1) w.r.t. +// 1/ instead of 1/. However, the mean number of photo-electrons can be +// expressed by with the relation: +// +// = c_phe * (2) +// +// Thus: +// +// 1/ = 1/c_phe * 1/ (3) +// +// and: +// +// Slope = F^2 / c_phe +// +// In the ideal case of having equal photo-multipliers and a perfectly flat-fielded camera, +// the "slope" -values should thus all be closely centered around F^2/c_phe. +// +TH2F *MCalibrationIntensityChargeCam::GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom) +{ + + TH2F *hist = new TH2F("hist","R vs. Inverse Charges - Fit results",45,-0.02,0.02,45,0.,30.); + hist->SetXTitle("Offset [FADC counts^{-1}]"); + hist->SetYTitle("F^{2} / / [FADC count / phe]"); + hist->SetFillColor(kRed+aidx); + + MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(); + + for (Int_t npix=0;npixGetSize();npix++) + { + + if (geom[npix].GetAidx() == aidx) + { + TGraph *gr = GetRazmikPlot(npix); + gr->Fit("pol1","Q"); + hist->Fill(gr->GetFunction("pol1")->GetParameter(0),gr->GetFunction("pol1")->GetParameter(1)); + } + } + return hist; +} + + Index: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.h (revision 5593) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.h (revision 5594) @@ -9,4 +9,12 @@ #include "MCalibrationChargeCam.h" #endif + +#ifndef MARS_MCalibrationCam +#include "MCalibrationCam.h" +#endif + +class TGraphErrors; +class TH2F; +class MGeomCam; class MCalibrationIntensityChargeCam : public MCalibrationIntensityCam @@ -36,4 +44,9 @@ ((MCalibrationChargeCam*)GetCam())->SetNumPhotonsPINDiodeMethodErr(f); } + TGraphErrors *GetRazmikPlot( const UInt_t pixid ); + TGraphErrors *GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col=MCalibrationCam::kNONE); + TGraphErrors *GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col=MCalibrationCam::kNONE); + TH2F *GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom ); + ClassDef(MCalibrationIntensityChargeCam, 1) // Container Intensity Charge Calibration Results Camera };