Changeset 5594


Ignore:
Timestamp:
12/14/04 18:18:01 (20 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r5593 r5594  
    2828  * msignal/MExtractTimeAndChargeSpline.cc
    2929    - fixed some smaller bugs for special cases
     30
     31  * mcalib/MCalibrationIntensityChargeCam.[h,cc]
     32    - implemented automatic plotting of Razmik plot and others
    3033
    3134
  • trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc

    r5046 r5594  
    4141#include "MCalibrationIntensityChargeCam.h"
    4242#include "MCalibrationChargeCam.h"
     43#include "MCalibrationChargePix.h"
     44
     45#include "MGeomCam.h"
     46#include "MGeomPix.h"
    4347
    4448#include <TOrdCollection.h>
     49#include <TGraphErrors.h>
     50#include <TH2F.h>
     51#include <TF1.h>
    4552
    4653ClassImp(MCalibrationIntensityChargeCam);
     
    7178}
    7279
     80// -------------------------------------------------------------------
     81//
     82// Returns a TGraphErrors with the number of photo-electrons vs.
     83// the extracted signal of pixel "pixid".
     84//
     85TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col)
     86{
     87 
     88  const Int_t size = GetSize();
     89 
     90  TArrayF phe(size);
     91  TArrayF pheerr(size);
     92  TArrayF sig(size);
     93  TArrayF sigerr(size);
     94 
     95  for (Int_t i=0;i<size;i++)
     96    {
     97      //
     98      // Get the calibration cam from the intensity cam
     99      //
     100      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
     101
     102      if (col != MCalibrationCam::kNONE)
     103        if (cam->GetPulserColor() != col)
     104          continue;
     105      //
     106      // Get the calibration pix from the calibration cam
     107      //
     108      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
     109      //
     110      // Don't use bad pixels
     111      //
     112      if (!pix.IsFFactorMethodValid())
     113        continue;
     114      //
     115      phe[i]    = pix.GetPheFFactorMethod();
     116      pheerr[i] = pix.GetPheFFactorMethodErr();
     117      //
     118      // For the calculation of Q, we have to use the
     119      // converted value!
     120      //
     121      sig   [i] = pix.GetConvertedMean();
     122      sigerr[i] = pix.GetConvertedMeanErr();
     123    }
     124 
     125  TGraphErrors *gr = new TGraphErrors(size,
     126                                     sig.GetArray(),phe.GetArray(),
     127                                     sigerr.GetArray(),pheerr.GetArray());
     128  gr->SetTitle(Form("%s%3i","Pixel ",pixid));
     129  gr->GetXaxis()->SetTitle("Q [FADC counts]");
     130  gr->GetYaxis()->SetTitle("photo-electrons [1]");     
     131  return gr;
     132}
     133
     134// -------------------------------------------------------------------
     135//
     136// Returns a TGraphErrors with the number of photo-electrons vs.
     137// the extracted signal over all pixels with area index "aidx".
     138//
     139// The points represent the means of the pixels values, while the error bars
     140// the sigma of the pixels values. 
     141//
     142TGraphErrors *MCalibrationIntensityChargeCam::GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col)
     143{
     144 
     145  const Int_t size = GetSize();
     146 
     147  TArrayF phe(size);
     148  TArrayF pheerr(size);
     149  TArrayF sig(size);
     150  TArrayF sigerr(size);
     151 
     152  for (Int_t i=0;i<size;i++)
     153    {
     154      //
     155      // Get the calibration cam from the intensity cam
     156      //
     157      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
     158
     159      if (col != MCalibrationCam::kNONE)
     160        if (cam->GetPulserColor() != col)
     161          continue;
     162
     163      //
     164      // Get the area calibration pix from the calibration cam
     165      //
     166      MCalibrationChargePix &pix = (MCalibrationChargePix&)(cam->GetAverageArea(aidx));
     167
     168      phe[i]    = pix.GetPheFFactorMethod();
     169      pheerr[i] = pix.GetPheFFactorMethodErr();
     170      //
     171      // For the calculation of Q, we have to use the
     172      // converted value!
     173      //
     174      sig   [i] = pix.GetConvertedMean();
     175      sigerr[i] = pix.GetConvertedMeanErr();
     176    }
     177 
     178  TGraphErrors *gr = new TGraphErrors(size,
     179                                      sig.GetArray(),phe.GetArray(),
     180                                      sigerr.GetArray(),pheerr.GetArray());
     181  gr->SetTitle(Form("%s%3i","Area Index ",aidx));
     182  gr->GetXaxis()->SetTitle("Q [FADC counts]");
     183  gr->GetYaxis()->SetTitle("photo-electrons [1]");     
     184  return gr;
     185}
     186
     187// -------------------------------------------------------------------
     188//
     189// Returns a TGraphErrors with the 'Razmik plot' of pixel "pixid".
     190// The Razmik plot shows the value of 'R' vs. 1/Q where:
     191//
     192//        sigma^2         F^2
     193//  R =   -------    =  ------
     194//         <Q>^2        <m_pe>
     195//
     196// and 1/Q is the inverse (mean) extracted signal
     197//
     198TGraphErrors *MCalibrationIntensityChargeCam::GetRazmikPlot( const UInt_t pixid )
     199{
     200 
     201  const Int_t size = GetSize();
     202 
     203  TArrayF r(size);
     204  TArrayF rerr(size);
     205  TArrayF oneoverq(size);
     206  TArrayF oneoverqerr(size);
     207 
     208  for (Int_t i=0;i<size;i++)
     209    {
     210      //
     211      // Get the calibration cam from the intensity cam
     212      //
     213      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
     214      //
     215      // Get the calibration pix from the calibration cam
     216      //
     217      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
     218      //
     219      // Don't use bad pixels
     220      //
     221      if (!pix.IsFFactorMethodValid())
     222        continue;
     223      //
     224      // For the calculation of R, use the un-converted values, like
     225      // in the calibration, since:
     226      //                 C^2*sigma^2     sigma^2
     227      //  R(lowgain) =   -----------  =  ------  = R
     228      //                 C^2*<Q>^2       <Q>^2
     229      //
     230      const Float_t mean = pix.GetMean();
     231      const Float_t meanerr = pix.GetMeanErr();
     232      const Float_t rsigma = pix.GetRSigma();
     233      const Float_t rsigmaerr = pix.GetRSigmaErr();
     234      r[i]    = rsigma*rsigma/mean/mean;
     235      const Float_t rrelvar = 4.*rsigmaerr*rsigmaerr/rsigma/rsigma + 4.*meanerr*meanerr/mean/mean;
     236      rerr[i] = rrelvar * r[i] * r[i];
     237      rerr[i] = rerr[i] <= 0 ? 0. : TMath::Sqrt(rerr[i]);
     238      //
     239      // For the calculation of 1/Q, we have to use the
     240      // converted value!
     241      //
     242      const Float_t q  = pix.GetConvertedMean();
     243      const Float_t qe = pix.GetConvertedMeanErr();
     244      oneoverq   [i] = 1./q;
     245      oneoverqerr[i] = qe / q / q;
     246    }
     247 
     248  TGraphErrors *gr = new TGraphErrors(size,
     249                                     oneoverq.GetArray(),r.GetArray(),
     250                                     oneoverqerr.GetArray(),rerr.GetArray());
     251  gr->SetTitle(Form("%s%3i","Pixel ",pixid));
     252  gr->GetXaxis()->SetTitle("1/Q [FADC counts^{-1}]");
     253  gr->GetYaxis()->SetTitle("\sigma_{red}^{2}/Q^{2}");     
     254  return gr;
     255}
     256
     257// -------------------------------------------------------------------
     258//
     259// Returns a 2-dimensional histogram with the fit results of the
     260// 'Razmik plot' for each pixel of area index "aidx" (see GetRazmikPlot())
     261//
     262// The results of the polynomial fit of grade 1 are:
     263//
     264// x-axis: Offset (Parameter 0 of the polynomial)
     265// y-axis: Slope  (Parameter 1 of the polynomial)
     266//
     267// The offset is a measure of how well-known the supposed additional contributions
     268// to the value "reduced sigma" are. Because a photo-multiplier is a linear instrument,
     269// the excess fluctuations are linear w.r.t. the signal amplitude and can be expressed by
     270// the proportionality constant F (the "F-Factor").
     271// Adding noise from outside (e.g. night sky background) modifies the recorded noise, but
     272// not the mean extracted signal, due to the AC-coupling. Thus, noise contributions from outside
     273// (e.g. calculating the pedestal RMS)have to be subtracted from the recorded signal fluctuations
     274// in order to retrieve the linearity relation:
     275//
     276// sigma(signal)^2 / mean(signal)^2  = sigma^2 / <Q>^2 = F^2 / <n_phe>               (1)
     277//
     278// Any systematic offset in the sigma(signal) will produce an offset in the "Razmik plot"),
     279// characterized by the Offset of the polynomial fit. Thus, in an ideal case, all pixels have their
     280// "offset" centered very closely around zero.
     281//
     282// The "slope" is the proportionality constant F^2, multiplied with the conversion factor
     283// phe's to mean signal (because the "Razmik plot" plots the left side of eq. (1) w.r.t.
     284// 1/<Q> instead of 1/<n_phe>. However, the mean number of photo-electrons <n_phe> can be
     285// expressed by <Q> with the relation:
     286//
     287//  <n_phe> = c_phe * <Q>                                                            (2)
     288//
     289// Thus:
     290//
     291// 1/<n_phe> = 1/c_phe * 1/<Q>                                                       (3)
     292//
     293// and:
     294//
     295// Slope = F^2 / c_phe
     296//
     297// In the ideal case of having equal photo-multipliers and a perfectly flat-fielded camera,
     298// the "slope" -values should thus all be closely centered around F^2/c_phe.
     299//
     300TH2F *MCalibrationIntensityChargeCam::GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom)
     301{
     302 
     303  TH2F *hist = new TH2F("hist","R vs. Inverse Charges - Fit results",45,-0.02,0.02,45,0.,30.);
     304  hist->SetXTitle("Offset [FADC counts^{-1}]");
     305  hist->SetYTitle("F^{2} / <n_phe>/<Q> [FADC count / phe]"); 
     306  hist->SetFillColor(kRed+aidx);
     307
     308  MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
     309
     310  for (Int_t npix=0;npix<cam->GetSize();npix++)
     311    {
     312
     313      if (geom[npix].GetAidx() == aidx)
     314        {
     315          TGraph *gr = GetRazmikPlot(npix);
     316          gr->Fit("pol1","Q");
     317          hist->Fill(gr->GetFunction("pol1")->GetParameter(0),gr->GetFunction("pol1")->GetParameter(1));
     318        }
     319    }
     320  return hist;
     321}
     322
     323
  • trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.h

    r5046 r5594  
    99#include "MCalibrationChargeCam.h"
    1010#endif
     11
     12#ifndef MARS_MCalibrationCam
     13#include "MCalibrationCam.h"
     14#endif
     15
     16class TGraphErrors;
     17class TH2F;
     18class MGeomCam;
    1119
    1220class MCalibrationIntensityChargeCam : public MCalibrationIntensityCam
     
    3644    ((MCalibrationChargeCam*)GetCam())->SetNumPhotonsPINDiodeMethodErr(f); }   
    3745 
     46  TGraphErrors *GetRazmikPlot( const UInt_t pixid );
     47  TGraphErrors *GetPheVsCharge( const UInt_t pixid, const MCalibrationCam::PulserColor_t col=MCalibrationCam::kNONE);
     48  TGraphErrors *GetPheVsChargePerArea( const Int_t aidx, const MCalibrationCam::PulserColor_t col=MCalibrationCam::kNONE);
     49  TH2F         *GetRazmikPlotResults( const Int_t aidx, const MGeomCam &geom );
     50 
    3851  ClassDef(MCalibrationIntensityChargeCam, 1) // Container Intensity Charge Calibration Results Camera
    3952};
Note: See TracChangeset for help on using the changeset viewer.