Changeset 5594 for trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc

Timestamp:

12/14/04 18:18:01 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc

@@ -41 +41 @@
 #include "MCalibrationIntensityChargeCam.h"
 #include "MCalibrationChargeCam.h"
+#include "MCalibrationChargePix.h"
+
+#include "MGeomCam.h"
+#include "MGeomPix.h"
 
 #include <TOrdCollection.h>
+#include <TGraphErrors.h>
+#include <TH2F.h>
+#include <TF1.h>
 
 ClassImp(MCalibrationIntensityChargeCam);
@@ -71 +78 @@
 }
 
+// -------------------------------------------------------------------
+//
+// 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;i<size;i++)
+    {
+      //
+      // Get the calibration cam from the intensity cam
+      //
+      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+
+      if (col != MCalibrationCam::kNONE)
+        if (cam->GetPulserColor() != 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;i<size;i++)
+    {
+      //
+      // Get the calibration cam from the intensity cam
+      //
+      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+
+      if (col != MCalibrationCam::kNONE)
+        if (cam->GetPulserColor() != 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 =   -------    =  ------
+//         <Q>^2        <m_pe>
+//
+// 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<size;i++)
+    {
+      //
+      // Get the calibration cam from the intensity cam
+      //
+      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+      //
+      // Get the calibration pix from the calibration cam
+      //
+      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
+      //
+      // Don't use bad pixels
+      //
+      if (!pix.IsFFactorMethodValid())
+        continue;
+      //
+      // For the calculation of R, use the un-converted values, like 
+      // in the calibration, since: 
+      //                 C^2*sigma^2     sigma^2
+      //  R(lowgain) =   -----------  =  ------  = R
+      //                 C^2*<Q>^2       <Q>^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 / <Q>^2 = F^2 / <n_phe>               (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/<Q> instead of 1/<n_phe>. However, the mean number of photo-electrons <n_phe> can be 
+// expressed by <Q> with the relation:
+//
+//  <n_phe> = c_phe * <Q>                                                            (2)
+//
+// Thus: 
+//
+// 1/<n_phe> = 1/c_phe * 1/<Q>                                                       (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} / <n_phe>/<Q> [FADC count / phe]");  
+  hist->SetFillColor(kRed+aidx);
+
+  MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
+
+  for (Int_t npix=0;npix<cam->GetSize();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;
+}
+
+